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EXECUTIVE  SUMMARY 


Interferometries  Inc.  has  investigated  the  phase  errors  on  single  antenna  NAVSPASUR 
data.  We  find  that  the  single  antenna  phase  errors  axe  well  modeled  as  a  function  of 
signal  strength  only.  The  phase  errors  associated  with  data  from  the  Kickapoo  transmitter 
are  larger  than  the  errors  from  the  low-power  trzinsmitters  (».e.,  Gila  River  and  Jordan 
Lake).  Further,  the  errors  in  the  phase  data  associated  with  the  Kickapoo  transmitter 
show  significant  variability  among  data  tziken  on  different  days. 

We  have  applied  a  quadratic  polynomial  fit  to  the  single  antenna  phases  to  derive 
the  Doppler  shift  zmd  chirp,  and  we  have  estimated  the  formal  errors  associated  with 
these  quantities.  These  formed  errors  have  been  parameterized  as  a  function  of  peak  signal 
strength  and  number  of  data  frames.  We  find  that  for  a  typical  satellite  observation  the 
derived  Doppler  shift  has  a  formal  error  of  -^0.2  Hz  and  the  derived  chirp  has  a  formal  error 
of  1  Hz/sec.  There  is  a  clear  systematic  bias  in  the  derived  chirp  for  targets  illuminated 
by  the  Kickapoo  transmitter.  Near-field  effects  probably  account  for  the  larger  phase  errors 
and  the  chirp  bias  of  the  Kickapoo  transmitter. 
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I.  INTRODUCTION 


In  an  earlier  study  entitled  The  Determination  of  NAVSPASUR  Phase  Difference 
Errors  (hereafter  Paper  I),  Interferometries  Inc.  modeled  the  errors  present  on  the  phase 
difference  between  antennas  at  a  given  receiver  site.  In  a  multi-element  interferometer 
array  such  as  that  used  at  NAVSPASUR  receiver  sites,  it  is  the  phase  difference  data 
which  allows  the  calculation  of  the  direction  cosines  to  the  target.  The  phase  difference 
error  model  which  resulted  from  our  previous  analysis  is  a  crucial  factor  in  estimating  how 
well  the  direction  cosines  and  cosine  rates  are  determined. 

In  addition  to  the  phase  difference  information,  there  is  also  information  content  in  the 
time  variability  of  the  phases  on  each  of  the  individual  antenna  elements.  From  this  data 
we  can  obtain  the  Doppler  shift  and  the  time  rate  of  change  of  the  Doppler  shift  (hereafter 
Doppler  and  chirp)  for  the  target  in  question.  In  the  present  work  Interferometries  Inc.  heis 
analyzed  NAVSPASUR  data  to  determine  the  single  antenna  phase  errors  and  examine  how 
well  these  data  can  be  used  to  determine  Doppler  and  chirp. 

A  second  Interferometries  Inc.  report  entitled  NAVSPASUR  System  Performance 
Analysis  (hereafter  Paper  II)  contained  an  examination  of  how  well  satellite  position  and 
velocity  could  be  determined  by  NAVSPASUR  eis  a  single  system.  It  is  demonstrated  in 
Paper  II  that  the  Doppler  measurement  plays  an  important  role  in  the  determination  of 
position  and  velocity.  Hence,  the  present  analysis  bears  directly  on  the  question  of  overall 
system  accuracy. 

The  research  reported  in  this  paper  was  done  to  complement  the  previous  work  and  to 
cover  topics  not  considered  in  those  papers.  There  were  three  major  goals  accomplished  in 
this  study.  First,  the  single  antenna  phase  errors  were  modeled.  Second,  the  formal  errors 
on  the  determination  of  Doppler  and  chirp  were  parameterized.  Finally,  the  data  were 
examined  to  determine  whether  systematic  effects  are  present  which  bias  the  determination 
of  Doppler  and/or  chirp. 

There  zu:e  two  major  sections  to  this  report.  The  first  section  presents  the  method 
of  analysis  and  the  second  section  contains  the  results.  There  are  three  sets  of  results 
presented  herein:  the  models  of  the  single  antenna  phase  errors,  a  parameterization  of  the 
uncertainties  in  Doppler  and  chirp,  and  an  assessment  of  the  accuracy  or  “correctness”  of 
the  values  of  Doppler  and  chirp. 

II.  ANALYSIS  METHOD 

In  this  section  we  present  a  discussion  of  the  analysis  method  used  in  this  study. 
We  first  discuss  the  processing  method  applied  to  each  observation  and  then  consider  the 
global  analysis  performed  on  each  of  three  data  sets. 

During  the  initial  data  processing  it  became  evident  that  there  were  significant  dif¬ 
ferences  in  the  pha.se  errors  between  data  associated  with  the  Kickapoo  transmitter  and 


1 


that  associated  with  the  two  low-power  transmitters.  In  contrast,  for  data  associated  with 
any  given  transmitter,  there  were  no  significeint  differences  between  the  phase  errors  at 
different  receiver  sites.  We  have  therefore  analyzed  the  Kickapoo  and  non-Kickapoo  data 
separately,  but  have  otherwise  combined  the  data  from  the  different  receiver  sites. 

The  analysis  of  the  data  sets  had  two  major  parts.  First,  the  single  antenna  phase  error 
model  was  determined.  Second,  the  uncertainties  in  Doppler  and  chirp  were  calculated 
along  with  the  Observed  minus  Calculated  (O-C)  values  for  Doppler  and  chirp. 

Three  sets  of  data,  supplied  by  NAVSPASUR  management,  were  processed  in  the 
course  of  this  study.  Table  1  provides  a  summary  of  these  data  sets.  The  total  number  of 
identified  observations  is  listed  in  the  Table,  together  with  a  breakdown  by  transmitter. 
Also  shown  is  the  number  of  observations  which  were  eliminated  at  each  stage  of  the 
processing.  Each  of  these  data  sets  contained  data  from  all  three  transmitting  stations 
and  either  four  or  five  receiving  stations. 

Only  observations  for  which  the  target  satellite  was  identified  were  processed,  since 
only  those  scans  had  the  predictions  which  are  required  to  calculate  O-C.  Further,  only 
observations  recorded  in  the  36  Hz  wide  “full  Doppler”  channel  were  analyzed  since  there 
were  not  enough  data  from  the  narrower  channels  to  provide  a  statistically  significant 
sample. 

Observation  Processing 

An  example  of  a  typical  satellite  observation  is  presented  in  Exhibit  1,  which  contains 
data  from  the  Hawkinsville  receiver  site  (receiver  station  #5).  The  first  four  lines  of 
Exhibit  1  contain  header  information  which  identifies  the  NAVSPASUR  catalog  number 
of  the  satellite,  the  Doppler  channel  in  which  the  data  were  recorded,  the  time  and  date 
of  the  observation,  the  number  of  phaise  and  amplitude  data  frames,  the  predicted  time  of 
fence  crossing  and  the  predicted  satellite  state  vector  based  on  the  NAVSPASUR  catalog, 
and  other  control  information.  Following  the  four  heruier  lines  there  are  up  to  55  data 
lines,  each  of  which  contains  an  average  power  level  and  the  measured  phases  for  each  of 
the  dipole  arrays  at  the  receiver  site.  The  average  power  levels  are  expressed  in  deciBels 
below  1  mWatt,  i.e.,  144  -144  dBm.  Therefore,  smaller  numbers  correspond  to  stronger 

signals.  This  average  power  is  commonly,  if  incorrectly,  referred  to  as  the  amplitude  and 
we  will  adopt  this  usage. 

For  receiving  stations  1,  3,  4,  and  6  there  are  twelve  dipole  arrays,  hereafter  antennas; 
stations  2  and  5  have  only  eleven  antennas.  The  phase  at  each  antenna  is  measured 
54.98  times  per  second  and  recorded  as  a  6-bit  integer.  The  6-bit  S2mpling  explains  why 
the  phases  in  Exhibit  1  are  integers  between  0  and  63.  The  phase  measurement  is  thus 
quantized  in  units  of  approximately  6  degrees.  Each  set  of  antenna  phases  and  associated 
amplitude  is  referred  to  as  a  data  frame. 

The  first  step  in  processing  each  observation  was  to  remove  the  phase  wraps  from 
the  data.  These  phaise  wraps  occur  since  the  phase  is  measured  modulo  2rr  radians.  The 
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observation  in  Exhibit  1  contains  at  least  four  such  phase  wraps  on  each  antenna.  In  order 
to  track  the  phase  as  a  function  of  time,  these  phase  wraps  must  be  removed. 

The  phase  wraps  are  removed  in  subroutine  DEROT  (see  Appendix  A).  The  data 
frame  with  the  strongest  signal  is  located  and  phase  cycles  are  removed  starting  from 
that  point.  The  slope  of  the  phase  curve  is  determined  using  the  frames  immediately 
preceding  and  following  the  peak.  The  slope  is  then  used  to  predict  the  next  phase  value 
and  rotation(s)  of  phase  axe  added  or  subtracted  so  that  the  corrected  phase  is  within  ±7r 
radians  of  the  prediction. 

If  the  slope  at  the  amplitude  peak  is  not  well  determined,  the  phases  are  not  corrected. 
DEROT  returns  an  error  flag  if  any  one  of  the  antennas  is  not  processed.  When  this 
happens,  the  entire  observation  is  discarded  and  the  next  observation  is  read  and  processed. 
As  can  be  seen  from  Table  1,  typically  15-25%  of  the  available  observations  are  flagged 
during  the  derotation  process  and  eliminated  from  further  analysis.  Almost  all  of  the 
discarded  scans  had  signal-to-noise  ratios  (as  defined  below)  of  10  dB  or  less. 

Data  from  18  August,  1988  has  been  processed  to  exeimine  how  many  of  the  sintennas 
in  a  flagged  observation  cire  typically  not  successfully  processed  in  DEROT.  For  Kickapoo 
observations,  approximately  77%  of  the  scans  had  4  or  fewer  “bad”  antennas  and  about  41% 
had  only  one  antenna  which  was  not  successfully  processed  in  DEROT.  The  corresponding 
statistics  for  the  non-Kick  data  were  85%  and  52%,  respectively. 

The  fjiilure  of  the  derotation  algorithm  appears  to  be  caused  in  large  part  by  the 
relatively  crude  phase  recovery  scheme  we  have  adopted.  No  attempt  has  been  made 
to  process  these  observations  by  using  a  more  sophisticated  phase  recovery  2ilgorithm. 
Due  to  the  large  qu2uitity  of  data  processed  here,  a  more  elegant  derotation  algorithm 
would  irapjose  a  prohibitive  computational  burden.  We  have  analyzed  the  phase  errors 
on  those  antennas  which  successfully  derotated  for  a  random  sample  of  flagged  scans  and 
are  convinced  that  we  have  not  biased  our  results  by  dropping  these  scans  from  further 
consideration. 

The  decision  to  discard  the  whole  scan  rather  than  continue  processing  using  only 
the  recovered  amtennas  was  made  to  reduce  bookkeeping  complications.  While  this  does 
eliminate  a  portion  of  the  data,  there  aire  sufficient  data  that  survived  all  of  the  processing. 
In  particular,  there  remain  sufficient  low  amplitude  data  in  our  sample. 

Each  antenna  is  processed  separately  in  a  single  call  to  DEROT  with  the  phase  units 
changed  to  radians.  The  array  of  phases  returned  from  DEROT  for  the  observation  of 
Exhibit  1  is  shown  in  Exhibit  2,  with  the  location  of  the  strongest  signal  indicated.  For 
the  frames  near  this  point,  the  line-to-line  change  in  phase  is  well  behaved  and  the  several 
antennas  show  roughly  similar  phase  curves. 

A  careful  examination  of  Exhibit  2  shows  an  additioneil  deficiency  in  the  DEROT 
subroutine.  For  the  first  antenna,  the  change  in  phase  between  frames  18  and  19  is  in 
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the  wrong  direction  (t.e.,  the  change  on  antenna  #1  is  positive  with  increeising  distance 
from  the  peak  data  frame,  while  the  other  10  have  negative  phase  changes).  The  DEROT 
cilgorithm  is  not  accurate  at  average  power  levels  of  ^  -153  dBm. 

Previous  work  (Paper  I)  has  indicated  that  the  effective  noise  level  for  “full  Doppler” 
mode  NAVSPASUR  observations  is  -152  dBm.  Therefore,  it  was  decided  that  -152  dBm 
would  be  used  as  an  amplitude  cutoff.  The  observations  were  edited  to  eliminate  those 
frames  with  lower  average  power.  This  editing  was  done  by  moving  forward  and  backward 
from  the  peak  amplitude  to  locate  the  first  frame  in  each  direction  with  power  <  -152 
dBm.  These  frames  and  all  data  frames  beyond  these  were  deleted. 

In  some  observations,  the  signal  falls  below  -152  dBm  as  one  moves  away  from  the 
pecik  and  then  comes  back  to  greater  strength.  Only  the  contiguous  frames  with  amplitude 
>  -152  dBm  surrounding  the  peak  were  used  in  further  processing.  Exhibit  3  contains  the 
observation  of  Exhibit  2  after  editing.  The  series  of  Exhibits  1-2-3  illustrates  the  above 
discussion.  In  Exhibit  1,  there  is  a  secondary  peak  at  frame  10.  However,  the  data  editing 
removes  the  first  twenty  frames. 

No  further  processing  was  done  on  an  observation  unless  10  or  more  frames  remained 
after  the  data  were  edited.  The  minimum  processing  requirement  of  ten  frcimes  was  ar¬ 
rived  at  by  trial-and-error.  It  was  determined  that  the  results  of  our  phase  fitting  were 
generally  not  reliable  for  those  observations  contauning  fewer  than  10  frames.  This  criterion 
eliminated  approximately  10%  of  the  Kickapoo  and  about  6%  of  the  non-Kick  observations. 

The  actual  pairameter  fitting  for  each  etntenna  was  done  in  subroutine  WGHTPOLY, 
which  is  a  generalized,  weighted,  least-squares  polynomial  fitting  routine.  There  are  two 
key  assumptions  in  this  procedure.  First,  the  phaise  errors  are  assumed  to  be  random  with 
a  Gaussian  distribution.  Second,  in  the  absence  of  noise,  the  data  and  the  model  would 
agree  exactly. 

Previous  results  have  indicated  that,  for  most  observations,  the  time  dependence  of 
the  phase  on  each  antenna  is  well  modeled  by 


(f>{t)  =  Co  +  Ci(t  —  <o)  +  ~ 


(1) 


where  is  the  phase  at  time  t,  and  cq,  ci,  C2  represent  the  fitting  coefficients.  The 
reference  time  <o  is  taken  to  be  the  predicted  pMs  time. 

The  constant  term  coefficient  cq  is  the  phase  at  the  reference  time  to-  At  t  =  to> 
Cl  =  d4)/dt  is  the  the  true  rate  of  change  of  the  phase.  In  units  of  rotations/sec,  this 
is  the  Doppler  shift  at  time  to  with  respect  to  the  center  frequency  of  the  digital  filter. 
The  Doppler  shift  with  respect  to  the  NAVSPASUR  transmit  frequency  is  obtained  by 
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summing  this  coefficient  with  the  center  frequency  of  the  digital  filter.  The  chirp  is  given 
by  |c2,  in  units  of  rotations/sec^. 

The  call  to  WGHTPOLY  requires  three  arrays:  time,  phase,  and  expected  phase  er¬ 
ror.  The  time  is  calculated  with  respect  to  the  predicted  pass  time  for  ease  in  comparing 
obser\’ations  and  predictions.  The  phase  is  the  DEROTated  data  for  each  individual  an¬ 
tenna.  The  expected  phase  error  is  predicted  from  the  amplitude,  and  must  be  known 
prior  to  running  the  program.  This  error  vector  controls  the  weighting  of  the  phase  data 
and  determines  the  uncertainties  in  the  parameters. 

The  fitting  routine  returns  the  three  polynomial  coefficients,  a  v'alue  of  the  reduced 
chi-square,  and  the  covariance  matrix.  The  reduced  chi-squau'e  returned  from  the  fitting 
routine  indicates  the  quadity  of  the  fit  for  the  parameter  determination.  It  is  calculated  by 

2  ^  [<^n  ~  (co  -r  Ci(tn  —  to)  -I- C2(tn  —  ^o)^)]^ 

^  El 

where  t„,  and  £„  are  the  values  of  the  phase,  time,  and  a  priori  error  for  the  data 
point;  N  is  the  total  number  of  data  points.  The  quantity  in  the  brackets  []  is  termed  the 
residual, 

A  value  of  reduced  chi-square  of  ~  1  tends  to  indicate  reasonable  agreement  between 
the  error  model  and  the  residuals  and  is  a  sign  of  good  agreement  between  the  data  euid 
the  polynomial  fit.  If  chi-square  is  small  {e.g.,  <  0.2),  it  may  indicate  that  the  errors  from 
the  error  model  are  significantly  larger  than  the  errors  on  the  data.  If  chi-square  is  large 
it  may  indicate  that  the  errors  from  the  error  model  are  too  smeill  and/or  the  data  is  not 
well  fit  by  a  second  order  polynomial.  The  various  possibilities  can  only  be  distinguished 
by  inspection  of  the  residuals. 

The  covariance  matrix  elements  (array  COVAR)  returned  from  the  fitting  routine 
allow  the  uncertaunties  (vector  UNC)  in  the  parameters  and  the  correlation  coefficients  to 
be  calculated.  The  parameter  uncertaunties  aire  given  by  the  square  roots  of  the  diagonal 
elements  of  COVAR,  and  the  correlation  coefficients  are  obtained  from  the  off-diagonal 
elements  scaled  by  the  uncertainties,  i.e.. 


UNC(n)  =  VCOVAR(n,n), 

CORL(n,m)  =  COVAR(n,m)/[UNC(n)  *  UNC(m)]. 


(за) 

(зб) 


The  fitting  routine  occaisionadly  returned  non-physical  values  of  the  parameters.  While 
it  is  not  clear  why  this  happened,  the  observations  with  “bad”  fits  were  easily  recognized. 
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These  solutions  had  calculated  Doppler  shifts  feir  outside  the  frequency  range  of  the  digital 
filters,  ve*y  large  values  for  the  reduced  chi-square,  and  large  correlation  coefficients.  The 
observations  were  eliminated  if  the  absolute  value  of  the  correlation  between  Cj  and  C3 
exceeded  0.99.  This  criterion  sucessfully  eliminated  all  the  bad  fits  without  eliminating 
any  useful  data.  Approximately  2%  of  the  Kickapoo  and  4%  of  non-Kick  observations  were 
eliminated  in  this  manner. 

The  parcimeter  uncertainties  £ire  formal  errors,  which  means  that  these  uncertainties 
represent  how  much  the  parameter  determination  will  be  affected  by  the  errors  in  the  data. 
The  value  of  the  uncertainty  represents  how  consistent  the  parameter  determinations  will 
be  from  one  antenna  to  another  at  a  single  site  (i.e.,  an  indication  of  the  precision  of  the 
measurement).  The  uncertainties  do  not  tell  whether  or  not  the  parameters  agree  with 
reality  (i.e.,  how  accurate  they  are). 

To  assess  the  accuracy  of  the  parameter  determination,  values  of  Doppler  and  chirp 
O-C  were  ctilculated.  NAVSPASUR  has  supplied  a  priori  positions  and  velocities  for  each 
satellite  observation  which  have  been  used  to  Calculate  the  expected  values  of  Doppler 
and  chirp.  These  expected  values  are  subtracted  from  the  Observed  values  derived  from 
the  parameter  fits,  to  obtain  O-Cs.  In  the  absence  of  systematic  errors  in  the  catalog 
predictions,  any  systematic  effects  which  bias  the  Doppler  and  chirp  determinations  can 
be  uncovered  by  analyzing  these  O-C  values. 

Error  Modeling 

The  uncertadnties  in  the  parameters  are  calculated  from  the  phase  error  model  only. 
While  the  goal  of  this  study  was  to  examine  how  well  the  Doppler  and  chirp  can  be 
meaisured,  a  necesseiry  first  step  was  to  determine  the  single  antenna  phase  error  model. 
A  proper  error  model  is  essential  to  properly  calculate  the  uncertainties. 

The  error  model  derived  in  Pap)€r  I  for  the  errors  on  the  phase  differences  was  used 
as  a  first-guess  error  model  in  this  analysis.  It  was  not  expected  that  the  single  antenna 
phase  errors  would  be  the  same  as  the  errors  on  differenced  phases  in  Paper  I.  First,  the 
errors  in  Paper  I  were  for  antenna  pairs  rather  than  for  individiial  antennas.  Second,  there 
were  systematic  effects  seen  in  Paper  I  which  led  to  roughly  consteint  phase  difference  error 
for  signal-to- noise  ratios  ^  lOdB.  Finally,  there  may  be  other  systematic  effects  present 
which  cancel  when  taking  phase  differences,  but  that  will  increase  the  errors  on  the  single 
antenna  phaises. 

A  version  of  the  fitting  program  called  ACCUM  was  written  to  accumulate  and  analyze 
the  phase  residuals.  The  residuals  are  the  phase  data  minus  the  vaJue  calculated  from  the 
fit  parameters.  In  the  absence  of  systematic  effects,  these  residuals  (the  term  inside  the 
brackets  [  ]  in  Equation  2)  will  be  due  to  random  measurement  errors  on  the  phases. 

In  ACCUM,  the  residuals  are  binned  by  amplitude  and  accumulated  for  all  successfully 
processed  data.  The  accumulated  residuals  are  used  to  calculate  the  mean  and  RMS  (Root- 
Mean-Square)  phase  error  at  each  amplitude.  Only  those  antenneis  for  which  the  value  of 
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the  reduced  chi-square  is  lower  than  a  preselected  value  are  used  in  order  to  eliminate 
“bad”  solutions.  The  program  also  totals  the  number  of  reduced  chi-square  values  <1.0 
amd  >  1.0. 

The  first  step  in  processing  each  data  tape  was  to  determine  the  phase  error  model. 
As  was  the  case  for  the  phase  difference  errors  in  Paper  I,  the  phase  error  statistics  here 
turn  out  to  be  well  described  by  a  single  parameter  model,  t.e.  they  depend  only  on  the 
received  aimplitude. 

The  correct  error  model  was  derived  by  analyzing  the  residuals  and  adjusting  the 
error  model  to  make  the  model  agree  with  the  RMS  of  the  post-fit  residuals.  An  initial 
guess  was  made  for  the  error  model  ?md  the  fitting  program  run.  The  RMS  of  the  post-fit 
residuals  was  used  to  revise  the  error  model  estimate,  and  the  entire  process  iterated. 

The  iteration  was  continued  until  two  criteria  were  met.  First,  the  number  of  an¬ 
tenna  solutions  with  reduced  chi-square  <1.0  was  approximately  equal  to  the  number 
with  reduced  chi-square  >  1.0.  Second,  the  calculated  RMS  error  at  each  amplitude  was 
consistent  with  the  predictions  of  the  error  model.  For  the  final  error  models  the  level  of 
agreement  was  usually  better  than  ±2®. 

Doppler  and  Chirp  Processing 

After  the  proper  phase  error  models  were  determined,  each  data  set  was  run  through 
the  program  FITS.  This  algorithm  re- processed  all  of  the  observations  and  generated  out¬ 
put  files  for  further  analysis. 

One  file  contained  the  results  of  the  paraimeter  determination  EUid  the  reduced  chi- 
square  for  each  antenna  of  an  observation  together  with  the  imcertainties  and  correlation 
coefficients  for  that  observation.  There  is  only  one  set  of  uncertainties  for  each  observation 
since  the  phase  error  model  depends  only  on  the  amplitude,  which  is  the  same  for  all 
antennas. 

The  second  major  output  file  contained  the  11  or  12  pairs  of  values  of  Doppler  and 
chirp  O-C  for  each  fully  processed  observation.  This  file  also  contains  the  reduced  chi- 
square  value  for  eiK:h  smteima  and  a  header  line  that  identifies  the  observation.  The  third 
major  output  file  lists  the  par2uneter  uncertainties,  the  i>eak  amplitude,  and  the  number 
of  contin  ■)us  data  frames  greater  than  -152  dBm. 

There  is  zdso  a  file  written  which  contains  a  line  for  eaw:h  observation  indicating  whether 
the  processing  succeeded.  The  FITS  algorithm  also  allows  a  file  to  be  written  which 
contains  all  of  the  post-fit  residuads.  This  has  not  been  done  on  any  of  the  complete  data 
sets  because  the  files  zae  too  large. 
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III.  RESULTS  and  ANALYSIS 

This  study  has  produced  three  sets  of  results:  the  phase  error  models,  the  parameter¬ 
ization  of  the  formal  Doppler  and  chirp  uncertainties,  and  the  analysis  of  the  Doppler  and 
chirp  residuals  (O-C). 

Phase  Error  Model 

The  data  tape  of  4Feb88  was  originally  processed  using  the  error  model  derived  in 
Paper  I.  It  was  clear  from  the  results  of  the  initial  processing  that  the  actual  phase  errors 
on  observations  involving  the  Kickapoo  transmitter  were  larger  than  both  the  input  error 
model  and  the  non-Kick  errors.  The  error  analysis  method  described  above  was  developed 
to  properly  determine  the  error  model  for  this  data  tape. 

There  were  two  clear  indications  that  the  initial  input  errors  were  too  small.  First, 
the  reduced  chi-squaxe  values  were  too  large.  Almost  all  of  the  chi-square  V2ilues  were  > 
1.0.  Second,  the  formal  uncertainties  in  the  parameters  were  much  less  than  the  observed 
variation  in  the  parameters.  In  particular,  the  chirp  determination  for  zdl  of  the  antennas 
should  be  the  same  to  within  the  formal  errors.  The  Kickap)oo  observations  had  antenna- 
to-antenna  variations  that  were  10-20  times  the  formal  errors. 

The  original  error  model  gave  more  reasonable  agreement  for  the  non-Kick  transmit¬ 
ters.  The  errors  (and  results)  for  Jordan  Lake  and  Gila  River  have  been  compared  with  no 
obvious  differences  seen  between  the  two  low  powered  transmitters.  While  there  has  not 
been  any  detailed  comparison  of  the  several  receiving  sites,  there  is  no  obvious  difference 
in  either  the  errors  or  other  results  for  the  different  receiving  stations. 

The  derived  error  models  for  the  three  data  sets  are  shown  in  Figures  la  and  lb. 
The  phase  error  in  degrees  is  plotted  versus  power  above  the  -152  dBm  noise  level.  The 
distribution  of  the  residuals  used  to  determine  the  error  model  has  been  examined  and  it 
is  in  reasonable  agreement  with  the  assumption  of  Gaussian  random  errors. 

The  Kickap>oo  errors  of  Figure  la  are  substantially  larger  than  the  non-Kick  errors 
of  Figure  lb.  Further,  the  Kickapoo  errors  for  the  three  tapes  are  not  in  agreement;  the 
three  curves  do  not  even  have  the  same  shape.  The  agreement  between  the  4Feb88  and 
25Apr88  tapes  at  large  signal-to-noise  (SNR)  may  not  be  significant  since  there  is  only  a 
small  amount  of  data  with  signals  this  strong  and  as  a  consequence  the  RMS  is  not  well 
determined.  There  were  several  iterations  done  for  each  Kickapoo  data  set  to  insure  that 
the  resulting  error  model  accurately  represented  the  errors  in  that  data  set.  The  error 
model  derived  from  one  tape  does  not  yield  an  acceptable  fit  to  the  other  data  tapes. 

All  three  sets  of  data  yield  essentially  identicad  error  models  for  the  non-Kick  transmit¬ 
ters  (see  Figure  lb).  Furthermore,  there  is  no  significant  difference  in  the  errors  between 
the  Jordan  Lake  and  Gila  River  transmitters. 


8 


The  errors  plotted  in  Figure  1  decrease  with  increasing  power,  eventually  levelling 
off  for  very  high  SNR  data  (^35  dB  SNR  for  Kickapoo  data,  ^20  dB  SNR  for  non-Kick 
data).  The  magnitude  of  the  errors  for  the  very  high  SNR  data  is  consistent  with  the  error 
expected  due  to  the  6  bit  phase  sampling,  which  imposes  a  roughly  6  degree  quantization 
on  the  reported  antenna  phase. 

It  was  expected  that  near-field  effects  would  lead  to  larger  errors  for  Kickapoo  obser¬ 
vations.  It  was  not  anticipated  that  the  difference  between  the  Kickapoo  and  non-Kick 
errors  would  be  so  dramatic.  Also,  there  was  no  reason  to  expect  the  marked  variation  we 
see  among  the  Kickapoo  error  models  for  the  three  different  data  sets.  The  scope  of  this 
study  has  not  allowed  a  quantitative  assessment  of  the  large  Kickapoo  phase  errors.  How¬ 
ever,  it  is  likely  that  the  larger  magnitude  of  the  Kickapoo  phzise  errors  is  due  to  near-  field 
effects.  Most,  if  not  aJl  satellites  observed  aje  at  ranges  shorter  than  the  far-field  distzmce 
of  s;  15,000  km.  In  the  near-field  region,  there  is  significant  structure  in  the  Kickapoo 
beam  which  may  account  for  the  larger  errors. 

Interferometries  Inc.  has  prepared  a  report  titled  Theoretical  Radiation  Patterns  of 
NAVSPASUR  Transmitter  Antennas,  which  demonstrates  the  dramatic  near-field  effects 
in  the  Kickap>oo  beam.  This  report  also  demonstrates  that  changes  in  the  bay-to-bay 
phasing  will  produce  laj*ge  changes  in  the  beam  pattern.  It  is  plausible  that  the  differences 
among  the  three  data  sets  are  caused  by  long-term  temporal  variations  in  the  bay-to-bay 
phasing  of  the  array. 

Doppler  and  Chirp  Uncertainties 

All  of  the  successfully  processed  observations  have  been  combined  so  that  the  uncer¬ 
tainties  in  the  Doppler  and  chirp  determinations  could  be  examined.  The  early  analysis 
made  it  clear  that  the  magnitude  of  the  uncert£unties  weis  determined  by  two  main  factors: 
the  peak  amplitude  and  the  number  of  data  frames. 

Tables  2a  and  2b  show  the  total  number  of  observations  for  the  Kickapoo  and  non- 
Kick  transmitters  as  a  fimction  of  number  of  data  frames  and  the  peak  amplitude.  It  was 
necessary  to  combine  the  three  data  sets  to  improve  the  quality  of  the  statistics. 

The  Kickapoo  observations  tend  to  have  larger  peak  signal  strength  and  more  data 
frames  thain  the  non-Kick  observations.  The  median  observation  from  Kickapoo  has  31 
data  frames  and  a  peak  signal  strength  of  -134  dBm.  The  corresponding  numbers  for 
non-Kick  are  28  data  frames  and  -141  dBm,  respectively.  (Note  that  the  maximum  of  55 
data  frames  per  observation  is  set  by  software  limitations  at  the  receivers.) 

Tables  3a-d  contain  the  meam  uncertaiinties  for  Doppler  and  chirp;  Tables  3a  amd  3b 
present  the  Kickapoo  results  while  Tables  3c  amd  3d  present  the  results  for  the  non-Kick 
tramsmitters.  Uncertaiinties  aire  only  reported  for  those  combinations  of  peak  amplitude 
and  number  of  data  frames  for  which  there  are  two  or  more  observations.  The  uncertainties 
decreaise  rapidly  as  the  number  of  data  frames  increases.  There  is  adso  a  significamt  decrease 
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in  the  uncertainties  as  the  peah  amplitude  increases,  with  some  flattening  out  for  very  high 
SNR. 

For  any  given  amplitude  and  number  of  data  frames,  the  uncertainties  in  Doppler  and 
chirp  are  significantly  smziller  for  non-Kick  observations  than  for  Kickapoo  observations 
(due  to  the  smaller  RMS  phase  errors).  However,  since  the  Kickapoo  transmitter  typically 
produces  more  data  frames  and  a  larger  peak  amplitude  than  non-Kick,  the  typiceil  obser¬ 
vation  will  have  roughly  equal  uncertainties  for  zuiy  of  the  three  transmitters.  The  median 
non-Kick  observation  has  uncertainties  of  0.14  Hz  in  Doppler  and  0.63  Hz/sec  in  chirp. 
The  corresponding  numbers  for  Kickapoo  are  0.24  Hz  emd  0.89  Hz/sec. 

Doppler  and  Chirp  Accuracy 

O-C  values  have  been  calculated  to  determine  the  degree  to  which  the  fitting  algorithm 
determines  the  “correct”  value  of  the  parameters.  Histograms  of  O-C  for  both  Doppler 
and  chirp  are  shown  in  Figures  2a-d.  The  histogram  for  the  Kickapoo  Doppler  (Figure  2a) 
is  more  sharply  peaked  than  the  corresponding  non-Kick  plot  (Figure  2b).  Both  Doppler 
histograms  peak  at  0  Hz,  which  indicates  that  there  axe  no  Izirge  systematic  offsets  in  the 
Doppler  determination.  If  the  distributions  are  taken  to  be  normal,  the  st2Lndard  deviations 
are  4.3  Hz  and  13.7  Hz  for  Kickapoo  and  non-Kick,  respectively. 

Mr.  F.  Lipp  of  NAVSPASUR  has  estimated  that  the  position  prediction  will  have  an 
RMS  error  of  —  8  kra.  In  Appendix  A  of  Paper  II,  the  partial  derivatives  of  Doppler 
and  chirp  with  respect  to  position  are  evaluated.  Typical  values  for  the  position  psu-tial 
derivatives  zire  ~  1  Hz/km.  This  implies  that  the  width  of  the  histogram  for  the  Kickapoo 
Doppler  may  be  due  to  the  quality  of  the  predictions.  A  more  accurate  assessment  of  the 
accuracy  of  the  observed  Doppler  will  require  the  analysis  of  NAVSPASUR  observations 
of  satellites  with  very  well  determined  a  priori  positions  and  velocities,  e.g.  LAGEOS, 
Starlette,  and  the  GPS  satellites. 

The  Doppler  O-Cs  for  the  non-Kick  observations  are  larger  than  those  for  Kickapoo, 
while  the  formal  uncertainties  are  smaller.  Furthermore,  the  non-Kick  histogram  (Figure 
2b)  is  nearly  flat  over  about^20  Hz;  the  distribution  is  clearly  not  Gaussian  and  is  much 
wider  thsui  predicted  by  the  formal  errors.  These  two  factors  may  indicate  the  presence  of 
one  or  more  small,  but  significant,  systematic  offsets  affecting  the  non-Kick  transmitter. 
Data  from  the  Jordan  L2dce  smd  Gila  River  transmitters  have  been  analyzed  separately, 
and  each  individual  station’s  data  closely  resembles  Figure  2b.  It  is  clear  that  further 
analysis  is  needed  to  verify  and  understand  these  possible  systematic  biases. 

Figures  2c  and  2d  present  the  histograms  of  chirp  O-C  for  the  Kickapoo  amd  non-Kick 
transmitters,  respectively.  The  chirp  O-Cs  for  the  non-Kick  observations  are  quite  smadl 
and  are  sharply  peaked  at  0  Hz/sec.  The  non-Kick  O-C  histogram  is  consistent  with  the 
uncertainties  of  Table  3d. 

The  situation  is  strikingly  different  for  the  Kickap>oo  chirp.  The  histogram  is  much 
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wider,  does  not  peak  at  zero,  and  is  strikingly  asyminetiic.  This  is  a  clear  indication  that 
there  are  systematic  effects  present  which  bias  the  Kickapoo  chirp  determination. 

This  bias  in  the  Kickapoo  chirp  has  been  previously  observed.  NRL  Memo  Re¬ 
port  #4831,  entitled  Expenmental  Observation  of  Naval  Space  Surveillance  Satellite  Sig¬ 
nals  with  an  Out-of-Plane  Receiving  Station,  by  S.  H.  Knowles,  W.  B.  Waltman,  and 
R.  H.  Smith  reported  an  average  chirp  bias  of  -1-10  Hz/sec  for  Kickapoo  observations. 

The  large  bias  on  the  Kickapoo  chirp  appears  to  be  due  primarily  to  near-field  effects. 
This  conclusion  has  been  tested  by  obtaining  additional  data  and  plotting  the  chirp  O- 
Cs  as  a  function  of  range.  Figures  3a  and  3b  display  the  chirp  O-C  versus  range  for 
Kickapoo  and  non-Kick  respectively  for  the  18  August,  1988  data  set.  The  magnitude  of 
the  Kickapoo  O-C  (Figure  3a)  is  clearly  a  function  of  range,  as  would  be  expected  for 
near-field  effects.  Figure  3b  shows  no  significant  range  effect.  It  is  clear  that  further  work 
is  necessary  to  quantitatively  understand  the  effect  of  range  on  the  chirp  determination. 

It  is  appropriate  at  this  point  to  discuss  how  well  we  expect  the  Doppler  shift  to  be 
determined  by  NAVSPASUR  in  normal  operation.  In  this  report,  we  have  parameterized 
the  formal  uncertainty  in  the  Doppler  for  a  single  antenna.  The  severed  antennas  at  a 
receiving  station  can  be  combined  to  reduce  the  formed  uncertainty  by  a  factor  of  about  3. 

It  is  clear  from  an  analysis  of  O-Cs  that  the  non-Kick  Doppler  observations  are  not 
as  accurate  as  the  corresponding  Kickapoo  observations.  Given  the  size  of  the  errors  in  the 
satellite  predicted  positions  we  have  been  unable  in  this  study  to  determine  the  absolute 
2w:c\iracy  of  either  the  Kickap>oo  or  non-Kick  Doppler  measurements.  Any  quantitative 
assessment  of  the  Doppler  accuracy  will  require  the  analysis  of  NAVSPASUR  observations 
of  satellites  with  very  well  determined  a  priori  position  and  velocity. 

It  is  our  opinion,  based  on  this  and  other  studies,  that  NAVSPASUR  should  be  able 
to  determine  the  Doppler  for  most  satellite  observations  to  an  accuracy  on  the  order  of 
0. 1-1.0  Hz.  Achieving  this  limit  will  require  that  the  systematic  effects  be  thoroughly 
understood.  F\irther  work  is  required  to  test  this  conclusion. 

rv.  CONCLUSIONS 

Interferometries  Inc.  has  investigated  the  phase  errors  on  single  antenna  NAVSPASUR 
data.  We  find  that  the  single  anteima  phaise  errors  are  well  modeled  as  a  function  of  signal 
strength  only.  The  phase  errors  on  observations  involving  the  Kickapoo  transmitter  are 
larger  than  the  errors  on  observations  from  the  low  power  transmitters.  FVirther,  the  three 
data  sets  which  were  anaJyzed  herein  yielded  signicantly  different  error  models  for  the 
Kickapoo  transmitter. 

The  formal  errors  in  the  determination  of  the  Doppler  shift  aind  chirp  have  been 
parameterized  as  a  function  of  peak  signal  strength  and  number  of  data  frames.  The 
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typical  satellite  observation  will  determine  the  Doppler  shift  with  a  formad  error  of  ~  0.2 
Hz  aind  the  chirp  with  a  formal  error  of  $  IHz/sec. 

The  Doppler  and  chirp  residuals  (0*C  values)  have  been  examined  for  indications  of 
possible  bias  due  to  systematic  errors.  The  distribution  of  both  the  Kickapoo  Doppler  and 
the  non-Kickapoo  chirp  residuals  axe  consistent  with  the  formal  errors  in  these  parameters 
and  the  RMS  errors  in  the  NAVSPASUR  position  and  velocity  predictions.  The  Kickapoo 
chirp  residuals  show  a  clear  ramge-dependent  bias  which  is  likely  to  be  caused  by  near¬ 
field  effects.  The  non-Kickapoo  Doppler  residuals  exhibit  a  non-Gaussian  distribution  and 
a  greater  RMS  deviation  than  expected,  given  the  formal  uncertainties  in  the  Doppler 
calculations  and  the  RMS  errors  in  the  NAVSPASUR  position  and  velocity  predictions. 
Further  study  is  needed  to  identify  the  source  of  these  larger  errors. 
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TABLE  1 

Summary  of  Observations 


Total  DEROT  <  10  Bad  Number  Fraction 

Observations  Error  Frames  Fit  Processed  Good 

4FEB88  13:30  to  4FEB88  17:00 

Kickapoo  5692  780  1133  200  3579  0.629 

non-Kick  4074  853  786  148  2347  0.576 


25MAR88  22:30  to  26MAR88  3:30 

Kickapoo  7191  1152  1081  321  4637  0.645 

non-Kick  5677  1219  700  208  3003  0.585 


21JUN88  15:30  to  21JUN88  22:00 

Kickapoo  8892  1214  1510  232  5856  0.659 

non-Kick  5677  1306  885  138  3348  0.590 


This  Table  lists  the  total  number  of  fuU-Doppler  observations  for  each  of  the  three  data 
tapes  processed  for  this  report.  The  observations  are  broken  down  between  the  Kickapoo 
and  non-Kick  transmitters  and  the  number  of  observations  eliminated  at  each  stage  of  the 
processing  is  listed.  Also  shown  are  the  total  number  of  observations  that  survived  the 
processing. 


THE  TOTAL  HUHBCK  ot  PEOCBSSED  OBSSEKVATIOIfS 
for  tho  EICKAPOO  TRARSMITTER 

••  a  PtmCTtOa  ot  PEAK  SICEAL  STAERGTN  and  tha  HUNBCR  OP  DATA  PRAHE  with  AMPLITUDE  >  -IS2  dBa 
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Table  2a  (Part  1) 
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TRI  TOTAL  MUNaiR  of  PROCESSCO  OSSSERVATIOMS 
foe  tho  KXCEAPOO  TRARSniTTCR 

••  •  PURCTIOa  of  PEAK  SXORAL  anfRaTK  aad  tho  RUMSER  OP  DATA  PRAME  with  AMPLITUDE  >  -152  dBa 


THE  TOTAL  HUHftCR  of  PROCESSED  OBSSERVATtOlfS 
for  tho  GXLA  RIVER  and  JORDAN  LAKE  TRANSMITTERS 

HAL  STRENGTH  and  tho  HUMBER  OP  DATA  FRAME  with  AMPLITUDE  >  -152 
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139  26  27  4  60  63  1  64  26  36  63  28 

141  18  20  0  46  66  61  49  18  29  69  25 

141  13  14  66  40  60  62  43  12  22  61  19 

142  2  6  49  33  44  47  36  3  17  47  16 

142  61  63  39  26  34  36  31  62  10  35  6 

144  60  50  30  17  28  30  19  62  0  28  63 

143  41  41  21  9  17  20  10  46  67  17  66 

146  30  30  13  61  6  6  0  30  43  6  47 

147  18  17  0  41  60  62  61  20  34  63  36 

147  12  6  62  36  46  47  46  11  26  46  25 


-5542.46 

-36.768 


EXHIBIT  1:  Raw  data.  The  first  three  lines  contain  header  ii^ormation  which  identifies  the  target 
satellite,  receiver  and  transmitter,  observation  time,  Doppler  channel,  predicted  pass 
time  and  state  vector,  and  some  control  information. 
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EXHIBIT  2:  The  phases  in  radians  as  output  by  DEROT.  The  locations  of  the  -153  dBm  amplitude 
cut-off  at  frame  number  20  (!)  and  the  amplitude  peak  at  frzune  number  39  (!!)  are 
shown.  Also  indicated  au’e  the  secondary  p>eak  at  frame  number  10  (*)  and  the  phase 
reversal  on  antenna  1  at  fr?ime  number  19  (♦•). 
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EXHIBIT  3:  Edited  data.  Only  the  data  frames  above  -153  dBm  surrounding  the  peah  stre  retained 
for  futher  processing.  A  careful  examination  of  the  phases  in  the  first  20  frames  in 
Exhibit  2  shows  that  this  is  a  reasonable  approach. 
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APPENDIX  A 


This  appendix  contains  annotated  listings  of  the  algorithms  used  in  this  study.  There 
are  two  versions  of  the  mmn  algorithm,  called  ACCUM  and  FITS.  Each  of  these  pro¬ 
grams  call  DEROT,  WGHTPOLY,  and  DPLCHR.  The  subroutine  GAUSSJ  is  czdled  by 
WGHTPOLY. 
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AUTHOR:  OR.  M.  ANDREWS 
DATE:  24  JUNE,  1988 

LANGUAGE:  FORTRAN  ANSI -77  (VAX/VMS  operating  system) 
FILE:  VX7770 : : SPACE : CANOREWS . CODE] ACCUM. FOR 


SUBROUTINES  CALLED:  OEROT.WGHTPOir 

LINK/LOAD  INSTRUCTIONS:  LINK  ACCUM,OEROT,UGHTPOLY. GAUSS J 

PROGRAM  DESCRIPTION:  This  routine  takes  out  the  phase  wraps  and  then 

does  a  weighted  third  order  polynomial  fit  to  the 
phase  as  a  function  of  time.  The  output  files  are 
set-up  to  test  the  error  model.  This  program  is 
almost  identical  to  FITS  except  for  error  analysis. 

PROGRAM  ALGORITHM  (PSEUDOCODE):  SEE  SUBROUTINES 

INPUTS  EXPLICIT:  There  are  three  input  files  required. 

1)  Error  model  table  which  is  transmit  station  specific. 

-  Read  from  FOR007. 

2)  Antenna  positions  in  Earth  centered  X,Y,Z  coordinate 

system.  Free  format  read  from  file  ANTPOS.DAT. 

•  Read  from  FOR008.  SPACE: UNDREWS. FILES] ANTPOS.DAT 

3)  NAVSPASUR  Phasei  data  files. 

-  Read  from  FQR009. 

IMPLICIT:  A  few  constants. 

OUTPUTS  EXPLICIT:  There  is  just  one  major  output  file.  F0R016  which 

contains  a  sunaary  of  the  errors  that  come  out 
of  the  fitting. 

MAJOR  VARIABLES:  In  order  of  appearance. 

ERR(45}  an  array  that  contains  the  model  of  phase  error  as  a 
fuKtion  of  amplitude.  The  erguaent  is  dB  above  -153  dBm. 

For  more  than  45  dB  use  the  value  at  45. 

CHICUT  the  CHI  SQUARED  cutoff.  Residuals  will  be  used  for 

sccumilation  in  the  error  analysis  only  if  the  value  of  CHI 
SQUARED  returned  by  WGHTPOlY  is  less  than  CHICUT. 

ISAT  NAVSPASUR  cstalogus  number  of  the  satellite. 

HRC  NTR  the  identifies  for  the  reciever  station  1-6  and  the 
transmit  station  7-9. 

IDOP  is  the  center  frequency  (  MINUS  10Hz  )  of  the  digital  filter 
w.r.t.  216.98  MHz  and  is  roughly  equal  to  the  Doppler  shift. 

lYR  IMON  IDAY  RSEC  are  the  year,  month,  day,  and  time  in  seconds 
on  that  day  of  the  start  of  data  record. 

NLINES  number  of  data  frames  in  that  observation. 

TPRED  is  the  predicted  pess  time  in  seconds  which  is  always 
close  to  RSEC. 

IAMP(55)  IPHS(55,12)  these  arreys  contain  tha  Phase4  data.  Thera 
will  be  NLINES  values  for  amplitude  and  NLINES  times  11  or 
12  phhses  (11  phases  if  NRC  ■  2  or  5,  12  otherwise).  The 
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units  on  lAMP  are  dB«  with  the  ainus  sign  OMITTED.  Smaller 
nunbers  are  stronger  si^ls.  IPHS  is  phase  sample  with 
6  bits;  hence  there  are  integers  fro*  0-63. 

RPHS(5S,12)  the  phases  returned  froai  DEROT.  The  rotations  of 
phase 

have  been  put  back  into  the  data.  The  units  are  radians. 

After  fitting,  this  array  is  changed  to  fold  the  residuals. 

lERR  the  error  flag  from  DEROT,  «  0  or  1.  If  IERR>1,  stop  processing 
of  this  observation  and  write  ■IERR>1*‘  to  »iit  13. 

MIAMP  the  maxinun  amplitude  in  the  observation,  e.i.,  the  mininun 
value  of  lAMP. 

IFRST  ILAST  NUM  these  variables  show  where  the  good  data  is  in  the 

observation.  IFRST  and  ILAST  identify  the  first  and  last  frames 
in  which  lAMP  <  1S3.  This  span  of  data  will  contain  the  peak 
amplitude.  NUN  is  the  nunber  of  good  frames.  If  NUM  <  10,  stop 
processing  this  observation  and  write  the  value  of  NUN  to  unit  13. 

AUSS)  A2(S5)  SIG(S5)  are  the  three  arrays  passed  to  UGHTPOLY.  They 

contain  the  times,  the  phase  data,  and  the  errors  on  the  phase. 

The  time  is  taken  w.r.t.  TPRED. 

A3<3)  contains  the  three  parameters  that  come  out  of  the  fit.  1  is 

a  phase  offset.  2  is  the  Doppler.  3  is  the  coefficient  of  T*T 
which  is  1/2  of  the  chirp. 

COVAR(3,3)  is  the  symmetric  covariance  matrix. 

CHISQ<12)  is  the  reduced  Chi  squared  returned  by  WCMTPOLY. 

N08S  INOFIT  NFEW  these  are  respectively  the  nudjer  of  observations 

in  the  data  set,  the  nuMber  for  which  DEROT  returned  IERR>1, 
and  the  nnbcr  in  which  HUM  <  10. 

RMEAN(66}  RMS<66)  are  the  mean  and  RMS  of  the  residuals  for  all  of 

the  accumjlated  data  at  a  given  amplitude. 

IACCUM(66,202}  this  array  is  used  to  accuaulate  the  residuals. 

MODIFIED:  19  July,  1988.  A  blank  read  added  due  to  change  in  the  content 
of  the  data  file. 

MODIFIED:  8  Septeaber,  1988.  Read  statements  changed  to  account  for  change  in 
data  foneat. 


IMPLICIT  REAL*8(A-H,0-Z) 

INTEGER*^  IAC0UM(66,202),HLES,NGRTR 
DIMENSION  X(99),T(99),A(20),RES(99),CHISa(12>, 

1  IAMP(S5),IPNS(S5,12),RPHS(5S,12>,ERR(AS), 

2  SIG(99),COVAR(3,3),IOCAN(66),RNS<66> 

TINCR  -  1./S4.98 
KNAX  -  0 
KMIN  ■  67 
INOFIT  ■  0 
N08S  ■  0 
NFEU  >  0 
NLES  ■  0 
NGTTR  ■  0 

Read  the  error  model  table. 

DO  J>1,A5 

REA0(7,*)  INTG,  ERR(J) 

ENOOO 
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Read  the  cutoff  value  of  Chi  Square  that  will  detenaine  which 
residuals  are  passed  to  lACCUM  Array. 

URITE(6.*)  'ENTER  CHI  SQUARE  CUTOFF' 

REAO(S.*)  CHICUT 

For  each  satellite  observation  read  the  data.  The  start  times 
is  TSTRT.  The  phases  get  loaded  into  the  X  array.  The  Y  array 
is  then  loaded  with  the  tisK  zero  at  first  occurance  of  maxiiius  AMPL. 
Call  OEROT  and  when  the  phases  successfully  derotate  call  UGNTPOLY 
NANT  times. 

301  REA0(9,*,EN0<321)  ISAT,NRC.NTR,IDOP,iyR,lN0N,10AY,RSEC, 

1  NLINES 

ISAT  «  ABS(ISAT) 

NOBS  -  NOBS  ♦  1 

Two  blank  reads  to  skip  lines  with  stuff  we  don't  need  and  TPREO. 

READ(9,*,END<321) 

REA0(9,*.EN0*321)TPRED 

REA0(9.*,EN0>321) 

NANT>12 

lF(NRC.Ea.2  .OR.  NRC.EO.S)  NANT>11 
00  1I«1, NLINES 

REA0(9.*.EN0*321)  lAMP(It),  (  tPHS(lI,J),  J>1,NANT) 

ENODO 

Find  line  nunber  of  first  occurance  of  maxinus  in  amplitude. 

NOTE:  largest  amplitude  is  smallest  ninber. 

MIAMP  «  190 

00  101  J«f. NLINES 

IF(  lANPfJI.GE.MIANP)  GOTO  101 
MIAMP  >  lAMP(J) 

NMAX  a  J 

101  CONTINUE 

Call  OEROT  which  will  derotate  the  phases  of  all  antemas  in  one  call. 

If  lERR  a  1  did  not  get  a  good  derotstion  so  go  to  next  observation. 

CALL  OEROT(NLINES,NANT,IANP,1PH$,RPNS,1ERR) 

IF(  lERR.EO.I)  INOFITalNOFIT  ♦  1 
IF(  lERR.EO.I)  GOTO  302 

Pass  only  the  "goo^  data  to  the  fitting  routine. 

Go  forward  and  backward  from  lina  NMAX  and  stop  if  lAMP  grastar 
than  152. 

IFRST  a  0 

DO  Ja1,NNAX-1 

IK  IANP(  NNAX-J)  .GT.  152)  GOTO  311 
ENOOO 
IFRST  a  1 

311  IF(  IFRST  .EG.  0)  IFRST  a  (MAX  •  J  «  1 

DO  JaNNAX^I. NLINES 

1F(  lAMP(J)  .CT.  152)  GOTO  312 
ENDOO 

312  lUST  a  J-l 

Oon't  continue  unless  there  are  10  or  more  good  lints. 

NUM  a  lUST  •  IFRST  *  1 
IF(NUH  .LT.  10)  NFEU  a  NFEW  *  1 
1F(NUM  .LT.  10)  GOTO  302 

Load  cp  the  SI6  array  by  using  the  ANDREWS  error  modal. 

Load  tp  X  with  times.  TaQ  at  first  occurance  of  max  asplituda. 


c 

DO  102  n«1,MUM 

J  •  H  ♦  IFRST  -1 
IMOX»  153  - 
IF(lNOX.GT.45)  1N0X-4S 
SIGdt)  >  ERRdNDX) 

Xdl)  >  (J  -  NMAX)*TINCR 

102  CONTINUE 
C 

C  Call  WGHTPOLY  for  each  antenna. 

C 

DO  103  II=1.NANT 
00  J>1.NUM 

JJ  »  J  ♦  IFRST  -  1 
Y(J)  >  RPHS(JJ,II) 

ENDOO 

C 

CALL  WGNTP0LY(X.Y.SlG.NUM.A,3.C0VAR,CHlSQdI>  > 

C 

IF(CHISO<n)  .LE.1.0)  NLES  «  NLES  *  1 
IF(CHlSCIdn  .'GT.1.0)  NGRTR  «  NGRTR  *  1 
C 

C  Calculate  the  residuals  and  replace  the  elenents  of  RPHS  so  that 
C  it  Mill  contain  residuals  rather  than  data. 

C 

DO  Jsl.NLINES 

T  *  (J  -  NNAX)*TINCR 
CALC  »  A(1)  ♦  A(2)*T  ♦  A(3)*T*T 
RPHS(J,I1)  =  RPHS(J.U)  '  CAl-C 
ENDOO 

103  CONTINUE 
C 

C  Increment  the  lACCUM  array  using  only  the  data  points  used  in  the  fitting, 

C  For  each  amplitude  bin  the  residuals  into  calls  .03  rsdians  wide. 

C  The  first  and  last  bins  contain  raspactivaly  those  residuals  lass  than 

C  -3.0  and  greater  than  «3.0  radians. 

C 

C  K  is  determined  by  the  magnitude  of  the  aay>lftude  and  INDEX  by  the 
C  magnitude  of  the  residual.  Amplitude  of  *152  dBm  is  the  smallest  and 
C  goes  into  the  K>1  bin.  Use  KNAX  to  keep  track  of  the  largest  as^lituda 
C  and  KMIN  for  smallest. 

C 

DO  104  L>IFRST,IUST 
K  «  lANPlL)  •  152 
IFIK.GT.O)  K>0 
K>  ASS(K)  *  1 
IFIK.GT.KMAX)  UMAX  -  K 
IFIK.LT.KNIN)  KMIN  ■  K 
DO  105  J-1,NAMT 

IF(  CNISQ(J)  .GT.  CHICUT  )  GOTO  105 

IF(  RPHSd.J)  .LT.  0.0  )  INDEX  ■  IHT((RPHS(L. J)-.015)/.03) 
IF<  RPHSd.J)  .GE.  0.0  )  INDEX  •  INT((RPHS(L. J>«.015)/.03} 
INDEX  •  101  *  INDEX 
IF(  INDEX  .LT.  1)  INDEX  >  1 
IF(  INDEX  .GT.  202)  INDEX  •  202 
IACCUM(K,INOEX)  ■  lACCUNd,  INDEX)  *  1 
105  CONTINUE 

104  CONTINUE 
C 

C  Return  to  reading. 

C 

302  REA0(9,*.EN0«321) 

GOTO  301 
C 

321  CONTINUE 

URITE(16,*)  'N08S«',N0BS.'  INOFIT>',INOFIT,'  NFEU-',NFEU 
URITE(16,*)  'N  Chisquara  LE  1>'.NLES,'N  Chisquara  GT  1.0«', NGRTR 
C 

C  When  the  end  of  the  file  is  reached  analyse  the  residuals  and  print 
C  out  the  results. 
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c 

C  Use  the  contents  of  tACCUM  to  calculate  the  MEAN  and  RMS  for  each 
C  aiiN)litude. 

C 

C  Print  the  results. 

C 

DO  106  KMIN.KMAX 
lAMPL  =11-153 
SUM  >  0.0 
SUM2  -  0.0 
NIN  =  0 

00  107  J=2,201 

IF(J  .LT.  101)  RESID  «  .03*(J-101>  -  .015 
1F(J  .GT.  101)  RESIO  «  .03*(J-101)  ♦  .015 
1F(J  .EQ.  101)  RESIO  «  0.0 
N  >  lACaJH(ll,J) 

NIN  =  NIN  4^  N 
SUM  =  SUM  4  N*RES10 
SUM2  «  SUM2  4  N*RES10*RES10 
107  CONTINUE 

IF(NIN.LT.IO)  GOTO  106 
RNEAN(ll)  =  SUM/NIN 

RMS(ll)  =  (  N1N/(N1N-1)  )*(  SUM2/N1N  -  RMEANCl 1 )*RMEAN(1 1)) 
RMS(ll)  *  SORT(  RMS(Il)  ) 

C 

UR1TE<16,201)  1AMPL,1ACCUM(11,1).1ACCUM(1 1,202), NIN, 

1  RMEAN(ll),RHS(It),  (  IACaJM(ll,K),  K=2,201) 

201  FORMAT! IX, 14,'  NL0W=',16,'  MM1CH«',16,'  N1N=' , 16, 10X, 

1  '  MEAN»<,F7.2,'  RMS=',F7.2,  20(/,10lS),/  ) 

106  CONTINUE 

ENO 
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PROGRM  FITS 


PROGRAM  FITS 


AUTHOR:  DR.  M.  ANDREWS 
DATE:  24  JUNE.  1988 

LANGUAGE:  FORTRAN  ANSI -77  (VAX/VMS  operating  tysten) 

FILE:  VX7770::SPACE:UNDREWS.COOE]FITS.FOR 

SUBROUTINES  CALLED:  DEROT. UGHTPOLY.DPLCHR 

LINK/LOAO  INSTRUCTIONS:  LINK  FITS.DEROT.WGHTPOLY.DPLCHR.GAUSSJ 

PROGRAM  DESCRIPTION:  This  routine  takes  out  the  phase  wraps  and  then  does 

a  weighted  third  order  polynonial  fit  to  the  phase  as 
function  of  time.  The  observed  Doppler  and  chirp  come 
out  of  this  fitting.  The  output  files  are  set-up  to  aid 
further  analysis. 

PROGRAM  ALGORITHM  (PSEUDOCODE):  SEE  SUBROUTINES 

INPUTS  EXPLICIT:  There  are  three  input  files  required. 

1)  Error  model  table  which  is  transmit  station  specific. 

-  Read  from  FOR007. 

2)  Antenna  positions  in  Earth  centered  X.Y.Z  coordinate 

system.  Free  format  read  from  file  ANTPOS.DAT, 

•  Read  from  F0R008.  SPACE: [ANDREWS. FILESIANTPOS. DAT 

3)  NAVSPASUR  Phase4  data  files. 

•  Read  from  F0R009. 

IMPLICIT:  A  few  constants. 

OUTPUTS  EXPLICIT:  There  are  five  output  files. 

11  NOTE  II  For  large  data  aets  it  is  necessary  to  restrict  the  output 
by  inserting  C,  comment,  at  the  appropriate  points. 

1)  F0R012  will  contain  a  header  line  to  identify  the 

observation.  The  derived  Doppler  and  chirp  along  with 
the  Chi  Squared  value  on  the  fit  for  each  antenna  is 
printed.  The  uncertainties  in  the  quantities  and  ,^e 
elements  of  the  covariance  matrix  are  printed.  ^ 

2)  F0R013  contains  a  header  line  which  identifies  the 

satellite  and  lists  the  pass  time.  For  each 
observation,  list  the  transmitter  and  receiver  nubers 
and  *0000  FIT"  if  everthing  worked.  Otherwise,  list 
reason  for  no  fit. 

3)  F0R014  lists  the  maximus  aeblituda  and  the  nimber  of 

data  lines  with  mqplitude  >  *153  dta  and  the 
uncertanties  in  Doppler  and  chirp. 

4}  FOR015  lists  the  0*C,  Observed  minus  Calculated  values 
for  Doppler  and  chirp.  Observod  values  from  the 
fitting,  calculated  values  from  DPLCHR. 

5)  F0R014  lists  the  residuals  from  the  potymonial  fit. 

"NOTE*  This  output  file  will  be  comparable  in  size  to 
data  file. 

MAJOR  VARIABLES:  In  order  of  appearance. 

ERR(45)  an  array  that  contains  the  model  of  phase  error  as  a  function  of 

amplitude.  The  arguemnt  is  dl  above  -153  dBm.  For  sure  than  45  dB 
use  the  value  at  45. 

RECPOS(6,12,3)  TRPOS(9,3)  these  two  arrays  contain  respectively  the 

antenna  positions  for  the  receivers  and  transmitters.  The  file 
SPACE:  [ANDREWS. FILESJANTPOS.DAT  has  been  set-tg)  so  that  the  free 
format  read  in  the  code  will  work.  The  arguements  of  RECPOS  are 
station  #  (1-6),  antenna  i  (1-12),  coordinate  index  (1-3:  1-X, 
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2-Y,3-Z.  The  trensfflit  station*  are  nuitered  7-9  as  identified 
by  NAVSPASUR.  The  second  index  is  coordinate  as  above. 

ISAT  NAVSPASUR  catalogue  nunber  of  the  satellite. 

NRC  NTR  the  identifies  for  the  reciever  station  1-6  and  the  transmit 
station  7-9. 

lOOP  is  the  frequency  of  the  digital  filter  w.r.t.  216.98  NHi 
and  is  roughly  equal  to  the  Doppler  shift. 

lYR  IHON  IDAY  RSEC  ere  the  year,  month,  day,  and  time  in  seconds  on 
that  day  of  the  start  of  data  record. 

NLINES  ncnber  of  data  frames  in  that  observation. 

TPRED  SX  SY  SZ  VX  VY  VZ  the  predicted  pass  tine  and  the  position  and 
velocity  of  the  satellite  at  tiaw  TPRED. 

TOBS  AEU  ANS  REU  RNS  FDOP  FRATE  the  solution  time  and  the  values  of 

the  of  the  EW  angle  and  the  NS  angle  in  radians,  the  rate  of  change 
of  sine  of  the  EU  and  NS  angles  in  radians/sec.,  the  fitted  Doppler 
and  chirp  as  determined  by  NAVSPASUR  at  time  TOBS. 

lAMPlSS)  IPHS<55,12)  these  arrays  contain  the  Phase4  data.  There 
will  be  NLINES  values  for  amplitude  and  NLINES  times  11  or  12 
phases  (11  phases  if  NRC  «  2  or  S,  12  otherwise).  The  inits  on 
lAMP  are  dBm  with  the  minus  sign  OMITTED.  Smaller  rsmbers  are 
stronger  signals.  IPHS  is  phase  saaple  with  6  bits;  hence  there 
are  integers  from  0-63. 

RPHS(55,12)  the  phases  returned  from  OEROT.  The  rotations  of  phase 
have  been  put  back  into  the  data.  The  units  are  radians. 

After  fitting,  this  array  fa  changed  to  fold  the  residuals. 

lERR  the  error  flag  from  OEROT,  •  0  or  1.  If  IERR>1,  stop  processing 
of  this  observation  and  write  "lERRal*  to  unit  13. 

MIAMP  the  maxiiSM  aaplitude  in  the  observation,  e.i.,  the  minissjs 
value  of  lAMP. 

IFRST  I  LAST  NUM  these  variables  show  where  the  good  data  is  in  the 
observation.  IFRST  and  ILAST  identify  the  first  and  last  frames 
in  which  lAMP  <  1S3.  This  span  of  data  will  contain  the  peak 
amplitude.  NUN  is  the  nmber  of  good  frames.  If  NUN  <  10,  stop 
processing  this  observation  and  write  the  value  of  NUN  to  unit  13. 

AKSS)  A2(5S)  SIG(SS)  are  the  three  arrays  passed  to  UGHTPOLY.  They 
contain  the  times,  the  phase  dsts,  and  the  errors  on  the  phase. 

The  tins  is  taken  w.r.t.  TPRED. 

A3(3)  contains  the  three  parmaeters  that  come  out  of  the  fit.  1  is 
a  phase  offset.  2  is  the  Doppler.  3  is  the  coefficient  of  T*T 
which  is  1/2  of  the  chirp. 

C0VAR(3,3)  is  the  symmetric  covariance  matrix. 

CHISQ  is  the  reduced  Chi  squared  for  this  particular  fit. 

DOP(12)  CHRP(12)  contain  the  Doppler  and  chirp  for  the  given  antenna 
nwber  as  calculated  in  OPLCNR. 

C0EFF(5,12)  holds  the  answers.  (1,J)  is  the  phase  offset  in  rotations 

at  T«0  for  antenna  nuaber  J;  (2,J)  is  the  Doppler;  (3,J)  is  the  chirp 
(i,J)  is  the  value  of  reduced  Chi  squared;  (5,1*3)  contain  the  formal 
errors  on  offset  Doppler  and  chirp;  (5,4-6)  contain  the  correlation 
for  offset-Doppler,  offset-chirp,  and  Doppler-chirp. 

0NC(12,2)  contains  the  Observed  Minus  Calculated,  0-C,  values  for  each 
antenna  of  Doppler  (J,1)  and  chirp  (J,2). 
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For  chirp,  0-C  is  simply  O.S*A3(3)  -  CHRP  with  irits  of  Hz/sec. 

For  Doppler,  0-C  *  (  IDOP  *■  10.0  ♦  A2(2)  )  -  OOP  with  inits  of  Hz. 

Note:  lOOP  *  10.0  is  the  center  frequency  of  the  digital  filter, 

NOBS  INOFIT  NFEU  these  are  respectively  the  niifber  of  observations 
in  the  data  set,  the  nunber  for  lAich  DEROT  returned  lERR^I, 
and  the  nuiber  in  which  HUM  <  10. 

HOOIFtED:  29  August  to  conform  to  a  new  data  format. 


IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  A1(55),A2(55),A3(3),COEfF<5,12), 

1  IAMP(55),IPHS(55,12),RPHS{55,12),ERR(45), 

2  SIG(55),COVAR(3,3),OOP(12),CHRP(12>, 

3  RECPOS(6.12.3),TRPOS(9,3),OMC(12,2) 

Initialize  constants. 

TINCR  *  1./54.98 
TUOPI  >  6.2831853 
INOFIT  X  0 
lOLD  >  0 
NOBS  >  0 
NFEU  >  0 
NFAIL  =>  0 

Read  the  phase  error  model  table. 

NOTE:  The  error  model  may  be  specific  to  a  tingle  transmitter 

and/or  receiving  channel.  The  data  must  be  presorted  to 
be  compatible  with  the  error  model. 

DO  J>1,45 

REA0<7,*)  into,  ERR(J) 

ENOOO 

Read  the  X,Y,Z  positions  for  each  receive  and  transmit  antenna. 

NOTE:  File  SPACE:  (AN0REUS.FaESlANTPOS.DAT  has  been  generated  so 
the  following  read  will  work. 

REA0(8,*)  RECPOS.TRPOS 

Read  the  data  and  print  a  headar  line  for  each  new  satellite. 

301  REA0<9,*,EN0>999}  ISAT,NRC,NTR.I00P,ITR,IN0N,IDAT,RSEC, 

1  NLINES 

IF(ISAT.EQ.O)  QOTO  302 
NOBS  ■  NOBS  ♦  1 
IFdSAT.EQ.lOlO)  GOTO  302 
URITE(13,201)  ISAT,irR,IN0N,l0AY,RSEC 
201  F0RNAT(1X,/,1X, 'SATELLITE  OBSERVED  ON  ',3I3,F10.3) 

lOLO  ■  ISAT 

302  CONTINUE 

Do  one  blank  read  and  then  read  position  and  velocity  and  the 
NAVSPASUR  solutions. 

NOTE:  This  is  the  correct  code  for  data  tape  of  ISAuguat. 

REA0(9,*, END-999) 

REA0(9,*,END>999)TPRE0,SX,SY,SZ.VX,VY,VZ 
READ(9, • , END-999} TOBS, AEW, ANS, REU, RNS, FDOP, FCHRP 

Read  NLINES  data  frames.  There  are  twelve  antennas  unless  the 
receiving  station  is  2  or  5. 
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NANT>12 

IFCMRC.EQ.Z  .OR.  NRC.EQ.5)  MANT*11 
DO  1I-1,NLINES 

READ(9.»,EH0=999)  IAMPCII),  (  IPHS(M,J),  J=1,NANT) 

ENDDO 

If  ISAT=0  the  observation  is  not  identified  and  should  not  be  processed. 
IfdSAT  .Efl.  0)  GOTO  303 

Find  line  nuifser  of  first  occurence  of  maxiiitJa  in  aitplitude. 

MOTE:  largest  a<Bplitude  is  smallest  nuiber . 


MIAHP  »  190 
DO  101  J=1,NLINES 

IF(  IA«P(J).GE.MIA#(P)  GOTO  101 
HIAMP  >  lAMP(J) 

NHAX  =  J 

101  CONTINUE 

Call  DEROT  Mhich  will  derotate  the  phases  of  all  antennas  in  one  call. 

If  lERR  *  1  did  not  get  a  good  derotation  so  go  to  next  observation  and 
print  a  line  out. 


CALL  0EROT(NLINES,NANT. lANP, IPHS,RPHS. lERR) 

IF(  IERR.E0.1)  THEN 
INOFITxINOFIT  ♦  1 

WRITE(13,202)  MRC.NTR.NLINES.MIAMP.RSEC 
202  F0RNAT(1X,4I5,F12.3,'  1ERR=1') 

GOTO  303 
ENOIF 


C 

C 

C 

C 

C 

C 


304 

C 


305 


Pass  only  the  “good-  data  to  the  fitting  routine. 

Go  forward  and  backward  from  line  MMX  end  stop  If  lAHP  greater 

Pwio^  results  indicate  *152  <Mm  is  the  effective  noise  floor. 


IFRST  •  0 

00  J»1,MMX*1 

IF<  IAM<  MMX-J)  .GT.  152)  GOTO  304 
ENOOO 
IFRST  •  1 

IF<  IFRST  .EO.  0)  IFRST  -  MMX  -  J  ♦  1 

00  J-MMX^I.NLIMES 

IF<  IAM<J)  .GT.  152)  GOTO  305 
ENOOO 

lUST  ■  J*1 


C 

C 

C 

C 

C 


203 


Don't  continue  unless  there  are  10  or  more  good  lines.  ^  ^ 
Previous  results  indicate  that  at  least  10  frames  are  needed 
for  reasonable  fits. 


1  •  lUST  -  IFRST  ♦  1 
IFCNUM  .LT.  10)  THEM 

URITE(13,203)  MRC,MTR,MlIMES,MIAM,RSEC,MUm 
F0RMT(1X,4I5,F12.3,'  MUM', ID 


MFEW  ■  MFEW  ♦  1 
GOTO  303 
ENOIF 


C 

C 

c 

C 

C 


Load  ip  the  SI6  srrsy  by  using  the  ANDREWS  error  model. 

Load  up  A1  with  times.  T>0  at  TPREO. 

TSTRT  *  RSEC  -  TPREO 
DO  102  II-1,NUN 

J  «  II  ♦  IFRST  -1 
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INDX*  153  -  lAMP(J) 

1F(INDX.GT.4S)  INOX-45 
SIGdl)  >  ERR(INOX) 

Al(II)  =  (J  -  0.5)*T1NCR  ♦  TSTRT 
102  CONTINUE 

C 

C  Call  UGHTPOLY  for  each  antenna. 

C  Call  DPLCHR  for  each  antenna. 

C  Results  with  proper  units  go  into  COEFF. 

C 

00  103  II-1.NANT 
DO  J-1,NUM 

JJ  -  J  +  IFRST  -  1 

A2  uiU  contain  the  data. 

A2(J)  -  RPHS(JJ,II> 

ENOOO 


CALL  UGHTPOLY(Al,A2,SIG,NUM,A3,3.COVAR,CHISa) 

The  matrix  inversion  will  fait  for  some  scans.  Test  the  value  of 
C0VAR(1,3)  to  see  if  the  result  is  reasonable. 

TEST  »  A8S(  C0VAR(1,3)/  SORT(  COVARtI, 1)»C0VAR(3,3)  )  ) 
IF{TEST.GT.0.99)  THEN 

URITE(13,204)  NRC.NTR.HLINES.HIANP,RSEC 
204  FORMATdX,  4I5,F12.3,'  INVERSION  FAILED') 

NFAIL-NFAIL  ♦  1 
GOTO  303 

ENOIF 

CALL  DPLCNR<SX,Sr,S2,VX,Vr,VZ,00P(II),CHtP(lI), 

1  RECPOS(NRC,IW1),ltECPOS(NRC,n.2),RECPOS(NRC,n.3), 

2  TRP0S(NTR,1).TRP0S(NTR.2),TRP0S(NTR,3)) 

Change  the  units  to  rotations  so  that  Dopier  will  be  Hertz  and  chirp 
Hertz/sec.  This  puts  TWPI  in  a  nuiter  of  places.  The  center  frequency 
of  the  digital  filter  is  lOOP  *  10.0. 

**N0TE:**  The  chirp  is  tMice  A3(3) 

This  is  obvious  and  follows  from 

Phi(t)  -  Phi<T»0)  ♦  Doppler*!  ♦  0.5*chirp*T*T 

COEFFd.II)  *  A3d)/TW0PI 
CQEFF<2,II)  -  A3(2}/TU0PI 

GMC(II,1)  -  (  IDQP  «  10.0  *  C0£FF(2.n>  >  '  OOP(II) 
COEFF(3,II>  ■  2.0*A3(3}/TW0PI 
0MC<II,2)  -  C0CFF(3,1I)  •  CHRPdl) 

C0CFF(4,II)  «  CHISQ 

Calculate  the  residuals  and  replace  the  eleaiants  of  RPNS  so  that  it 
will  contain  the  residuals.  Residuals  will  be  in  units  of  radians. 

00  J«IFRST,ILAST 

T»  (J  -  0.5  )*TINCR  ♦  TSTRT 
CALC  •  A3<1)  ♦  A3<2)*T  ♦  A3<3)*T*T 
RPNS(J,M)  ■  RPNS(J,It>  -  CALC 
ENOOO 

103  CONTINUE 

Load  the  isKertainties  and  the  correllation  coefficients  into 
C0EFF(5,1-6). 


C0EFF(5,1)  >  SQRT(  COVAR<1,1)  )/TW0PI 
CQEFF(5,2)  >  S0RT(  C0VAR(2,2)  )/TW0PI 
C0EFF(5,3)  •  2.0*SORT(  COVAR(3,3)  )/TM0f>I 
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COEFF(5,4)  =  C0VAR(1,2)/  SO»T(  COVAJJ{1,1)*COVAR(2,2)  ) 

COeFF(5,5)  »  C0VAR(1.3)/  SORT(  COVARd ,  1  )*COVAR(3,3)  ) 

COeFF(5,6)  =  C0VAR<2.3)/  SQRT{  COVAR(2,2)*COVAR(3.3)  ) 

C 

C  Print  out  the  results. 

C  Header  line  and  coefficients  go  into  12. 

C  File  13  will  contain  a  header  line  and  indicate  a  good  fit. 

C  File  U  has  uncertainties  in  Doppler  and  Chirp. 

C  File  15  0-C  for  Doppler  and  Chirp. 

C  File  16  the  residuals  front  the  fitting. 

C 

URITE(12,205)ISAT,NRC,NTR,1YR,IMON,IDAY,TPRED,NUM,MIAMP,IDO(>,RSEC 

205  FORHAT(1X,/,1X,17,5I3,F10.3,2I4.I7.F10.3) 

C 

C  For  reduced  output  add  connents  here. 

C 

DO  JJ=1,NANT 

UllITE(12.206)JJ,COEFF(1.Ja),COeFF(2.JJ),COeFF(3,JJ),CO£FF{4,JJ) 

ENDDO 

206  FOfiMAT(1X,I6,F8.2,  3F12.2) 

C 

WRITE{12,207)  <C06FF(5,K),  K*1.6) 

207  F0«MAT{1X,6F10.4) 

C 

URtrE(13,208}  NRC.NTR.MIAMP.NUM.RSEC 

208  FO«MAT<1X,415,F12.3,'  GOOD  FIT') 

C 

UI)ITE(14,209)  NRC,NTR,MtAMP,NUM.CaEFF(5,2),COCFF(5,3) 

209  FORMAT(1X,4I4,2F10.2) 

C 

C  For  reduced  output  add  conmenta  here. 

C 

»«;T6<i5,210)ISAT,  WC.tITR,  lYR,  JJOi,  I0Ay,TPRED,RSEC,HIAI»P,MUN, 

1  IDaf>,Daf>(NANT),CMRP(HANT) 

DO  L«1,MANT 

UIIITE{15,211)L,0MC{L,1),0NC(L,2>,C0CFF(4,L) 

ENOOO 

210  F0SNAT<1X,/,1X,I7,2I3,3X,3I3,/,2F12.3,4X,2I4,3X,I6,2X,2F8.2) 

211  F0RNAT(1X,I2,3F10.2) 

Print  out  header  tine  and  the  reeiduata. 

This  will  be  a  large  output  file  *  edd  conwents  to  reduce  outputi 

UltITE(16,212)  ISAT,NRC.NTR,TPRE0,RSEC 
DO  JJ-IFRST.ILAST 

URITE(16,213)  lAMPtJJ),  (  RPHS(JJ.K),K<1,NAIIT  ) 

EHDOO 

212  F0»MT(1X,/,1X,3I5.2F10.3) 

213  FORNAT(1X,I6,12F7.2) 

Return  to  reading. 

303  READ<9,*,END>999) 

GOTO  301 

Branch  to  this  point  when  EOF  encountered  during  read. 

999  CONTINUE 

M1TE(13,*)  'N08S»',N0iS,'  INOFIT»',INOFIT,'  NFEV'.NFEW, 

1  '  HFAIL«',HFAIL 

C 

END 
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SUBROUTINE  OEROT(NAMPL,NANT, lAMPL, tPHS.RPNS, lERR) 


rM* 


•••*«********************t 

>  SUBROUTINE  DEROT 
r 

A******************************************************************** 


AUTHOR:  DR.  H.  ANDREWS 
DATE:  24  JUNE,  1988 

LANGUAGE:  FORTRAN  ANSI -77  (VAX/VMS  Operating  system) 

FILE:  VX7770 : : SPACE : [ANDREWS . CODE] DEROT . FOR 

CALLING  ROUTINE:  FITS 

PROGRAM  DESCRIPTION:  This  a  somewhat  simple  minded  program  whcih  removes 

the  phase  wraps  from  NAVSPASUR  Phase4  data.  The 
first  occurance  of  the  maximua  amplitude  is  located, 
the  slope  of  the  phase  curve  is  determined,  and  the 
slope  used  to  predict  and  correct  the  phase. 

INPUTS  EXPLICIT:  All  contained  in  the  calling  list. 

OUTPUTS  EXPLICIT:  All  contained  in  the  calling  list. 

MAJOR  VARIABLES: 

NAMPL,  NANT  the  nunber  of  data  frames  and  the  nuiber  of  antennas. 

lAMPLlSS),  IPHS(SS,12)  arrays  that  contain  the  raw  integer  amplitudes 
and  phases. 

RPH$(5S,12)  the  real  nunter  values  of  the  phases  returned  by  this 
program.  Units  are  radians. 

lERR  the  error  flag  which  indicates  whether  or  not  all  the  phases 
have  been  successfully  processed.  ■  1  indicates  error. 


MODIFIED: 


IMPLICIT  REAL'S  (A-H,0-Z) 

DIMENSION  IAMPL(55},IPHS(55,12).RPHS(55,12) 

Set  constant  values.  FACTR  converts  from  six  bit  integers  to  radians. 

TW0PI«6.2831855 
FACTRw  TWOPI/M. 

IERR-0 

Reset  the  RPHS  array  to  0. 

DO  11-1,55 
DO  JJ-1,12 

RPHS(II,JJ)-0. 

ENOOO 
ENODO 

Find  the  index  of  the  minimua  of  the  amplitude  i.e.,  STONGEST 
signal.  Will  derotate  phases  by  starting  there  and  working 
out. 

ININ  ■  200 

DO  111  LL-1, NAMPL 

IFdAMPLlLL)  .GE.  ININ)  GOTO  111 
ININ  ■  lAMPL(LL) 

IREF  >  LL 
111  CONTINUE 

C 
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C  If  the  first  or  last  amplitude  has  been  selected  move  in  one  place. 

C 

IFfIREF.EO.I)  IREF*I*EF  ♦  1 
IF(IREF.EQ.NAMPL)  tREF=IREF  -  1 

Select  antenna  nunber  III. 

Unwind  the  phase  by  starting  at  IREF  index  and  going  to  top  and 
bottom. 

At  the  start  use  the  points  before  snd  after  IREF  to  calculate  a 
slope.  Use  the  slope  to  predict  the  next  value.  Adjust  the 
rotations  as  needed  to  get  close  to  the  predictor. 

Convert  phases  to  radians. 

DO  101  III-1,NANT 

IDIF1  =  IPHS(IREF,II1)  -  IPHS{IREF-1,ni) 
lF(IDIF1.LT.-32)  I0IF1  *  lOIFl  ♦  64 

IFdOIFI.GT.  32)  IDIF1  *  I0IF1  -  64 

IDIF2  »  IPHS(IREF*1,III)  -  IPMSdREF.IIl) 

IFdDIF2.LT.-32)  I0IF2  «  I0IF2  ♦  64 

IFdDIF2.CT.  32)  I01F2  »  ID1F2  -  64 

Check  to  see  that  these  differences  are  either  small  or  of  the  same 
sign. 

If  the  slope  is  not  well  determined  at  the  amplitude  peak,  set  lERR^I 
and  do  not  attempt  to  process  the  phases. 

IF(  ABSdOIFI  «  IDIF1)  .LT.  10)  GOTO  301 

IF(  ABSd0IF1)>ABSd0IF2)  .NE.  IABS(  I0IFUIDIF2  )  )  GOTO  961 

301  CONTINUE 

Set  ISLOPE  to  the  slope  at  peak  and  convert  phase  at  peak  to  radian*. 

ISIOP  ■  (  IDIF1  «  I0IF2  )/2 
RPHSdREF.III)  «  IPHSdREF,ttI)*FACTR 
NROT  >  0 
ISLOPE-ISLOP 

Work  from  peak  to  the  beginning  of  the  scan. 

DO  JJ«1,IREF-1 
1I>IREF-JJ 

Predict  the  phase  based  on  the  previous  phase  value  snd  the  slope. 
IPROCT  ■  IPHSdl^l.IM)  -  ISLOPE 

Calculate  difference  between  observed  and  predicted  integer  phase  value. 
Add  or  subtract  a  rotation  to  make  the  data  closer  to  the  prediction. 

lOIFF  ■  IPHSdl.llI)  -  IPROCT  IFdOIFF.LT.-32) 

NROT  •  NROT^I 

IFdOIFF.GT.  32)  NROT  >  NROT-1 
Calculate  phase  in  radians. 

RPHSdI.IlI)  »  (  IPHSdl.III)  ♦  0.5)*FACTR  ♦  NROT*TWOPI 


Recalculate  the  slop*. 

ISLOPE  «N1NT(  (  RPHSdU1,III)  -  RPHSdI,III)  )/FACTR) 
ENOOO 


SO 
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Do  the  same  thing  from  the  peak  to  the  end  of  the  scan 
Reinitialize  NROT  and  ISLOPE. 


ISLOPE=ISLOP 
NROT  «  0 

DO  II::1REF«1,NAMPL 

IPROCT  s  IPHS(I1-1,III)  ♦  ISLOPE 
IDIFF  >  1PHS( 11,111)  -  IPRDCT 

IF(I01FF.Lr.-32)  NROT  «  NROT+1 
IFdDIFF.GT.  32)  NROT  «  NROT-1 

RPHSdl.IM)  *  <  IPHSdl.IlI)  ♦  0.5)*F«TR  ♦  NROT*TMOPI 

ISLOPE  >  NINT(  (  RPHSdI.III)  -  RPHSd  1-1 , 1 1 1 )  )/FACTR) 
ENOOO 

101  CONTINUE 

RETURN 

981  tERR>1 

If  this  point  is  reached  the  slope  is  rwt  properly  determined  and 
no  predictor  is  calculated. 

RETURN 

END 


51 


c%  f>  t>  oooo  ooooooooooorioooooooooooooooooooorioooooooooorjoooooooo 


SUBROUTINE  UGHTPOLY(X,Y,SIG,NOATA.A.NA.COVAR,CHISQ) 


SUBROUTINE  WGHTPOLY 


AUTHOR:  DR.  M.  ANDREWS 
DATE:  24  JUNE,  1988 

LANGUAGE:  FORTRAN  ANSI -77  (VAX/VNS  operating  system) 

FILE:  VX7770: :SPACE: [ANDREWS. COOEIWGHTPOLY. FOR 

CALLING  ROUTINES:  ACCUM.FITS 

SUBROUTINES  CALLED:  GAUSSJ 

PROGRAM  DESCRIPTION:  This  is  a  general  weighted  polynomial  fitti^ 

routine.  The  fivwtional  form  it 
Y  »  A(0)  ♦  A(1)  ♦  A(2>*X**2  ♦  etc. 

The  routine  uses  Chi  Squared  minimization  to  determine 
the  polynomial  coefficients.  The  covariance  matrix 
and  the  reduced  Chi  Squared  are  returned. 

This  routine  is  adapted  frcm  an  algorithm  in  Nunerical 
Recipes  which  should  be  consulted  for  details. 

INPUTS  EXPLICIT:  All  in  calling  list 

OUTPUTS  EXPLICIT:  All  in  calling  list 

MAJOR  VARIABLES: 

X,Y,S  are  three  arrays  of  length  NDATA  which  contain  the  variable 
and  the  error  values. 

NDATA  is  the  number  of  data  points. 

A  is  an  array  of  length  NA  that  contains  the  fit  parmaeters. 

NA  is  the  nudser  of  parameters  solved  for.  The  polynomial  order 
is  NA-1. 

COVAR  is  two  diiasnsionsl  array  of  dimension  NA*NA.  The  diagonal 
elements  are  the  square  of  the  UKsrtainties  in  the  parameters 
The  off-diagonal  elements  are  the  correlation  coefficients 
scaled  multiplied  by  the  diagonal  elements. 

CHISQ  is  the  value  of  the  reduced  Chi  Squared. 


MODIFIED: 


IMPLICIT  REAL*8  (A-H,0-Z) 

The  dimension  of  the  arrays  is  the  calling  list  is  controlled  by  the  values 
of  NA  and  NDATA.  Dimension  on  FX  and  S  limit  order  to  20  i.e.  X**19. 

DIMENSION  X(NDATA),Y(N0ATA),SIG(N0ATA).FX(20),S(20), 

1  A(NA),COVAR(NA,NA) 

Initialize  the  necessary  matrices. 

DO  J>1,HA 
DO  K«1,NA 

COVAR(J,K)  -  0.0 
ENDOO 

A(J)  *  0.0 
ENDOO 
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Loop  over  the  data  and  increment  the  normal  equations. 

NOTE:  Replace  a  zero  value  of  X  by  1*10-5  to  prevent  errors. 

DO  101  n>1,NDAT* 

DATA  «  X(II) 

1F(  ABS(  X(II)  )  .LT.  0.00001)  DATA  -  0.00001 
Load  the  array  FX(I)  by  FX(I)  *  0ATA**(I-1), 

DO  J>1,NA 

FX(J)  *  DATA**<J-1) 

ENDOO 

YM  =  Y(ll) 

SIG2I  «  1./SIG<II)**2 
DO  101  J^I.NA 
UT  «  FX(J)*SIG21 
DO  K>1.J 

COVARIJ.K)  «  COVAR(J,K)  «  WT*FX<K) 

ENDDO 

A(J)  «  A(J)  ♦  YM*WT 

101  CONTINUE 

Fill  in  across  the  diagonal  from  synmetry. 

IF(NA.Ea.l)  GOTO  102 
DO  102  J>2.NA 
DO  102  K«1,J-1 

C0VAR<K,J)  «  COVAR<J,K) 

102  CONTINUE 


Rescale  the  matrix  to  insure  nuaerical  stability. 

00  J*1,NA 

S(J)  «  SQRTl  C0VAR(J,J)  ) 

IF(  $(J)  .EQ.  0.0)  S(J)  •  1.0 
ENOOO 
DO  J>1,NA 

A<J)  ■  A(J)/S(J) 

DO  K<1,NA 

CCLARfJ.K)  ■  COVAR(J,K)/(  S(J)*S(K)  ) 

ENOOO 

ENDOO 

Call  GAUSSJ  to  invert  the  eaitrix. 

CALL  GAUSSJ(COVAR,NA,NA,A.1,1) 

Scale  back  to  correct  units. 

DO  J>1,NA 

A<J)  -  A(J)/S(J) 

DO  Kal.HA 

COVAR<J,K)  >  COVAX(J,K)/(  S(J)*S(K)  > 

ENOOO 

ENOOO 

Calculate  chi-squaro  by  suRtaino  the  differences  of  data,  Y's, 
and  the  values  calculated  using  the  fit  parameters. 

Convert  to  reduced  Chi  Squared  by  dividing  by  NDATA  -  NA  which 
is  the  nuiter  of  degrees  of  freedom. 

CNlSQaO.O 
DO  103  1I-1,N0ATA 
DATA  >  Xdl) 

IF(  ABS(  X(1I)  )  .LT.  0.00001)  DATA  ■  0.00001 
SUM«0.0 

DO  J>1,NA 

SUM  •  SUM  ♦  A{J)*0ATA**(4-1) 

ENOOO 
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RESID  «  Ydl)  -  sm 

CHISQ  «  CHISQ  ♦  (  RESlO/SIOdI)  >**2 


103 


CORTIHOE 

CHISQ  »  CHISO/(  MOATA  -  MA  ) 

RETURN 

END 
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SUBROUTINE  I>PLCHR(X,V.2,VX,VY,V2,00PLR,CNIRP,XR,YR,ZR,XT,rT,ZT) 


SUBROUTINE  DPLCHR 


AUTHOR:  DR.  H.  ANDREWS 

DATE;  24  JUNE,  1988 

LANGUAGE:  FORTRAN  ANSI -77  (VAX/VMS  operating  system) 

FILE:  VX7770 : : SPACE : [ANDREWS . CODE] DPLCHR . FOR 

CALLING  ROUTINE:  FITS 

PROGRAM  DESCRIPTION:  This  subroutine  uses  the  position  and  velocity  of 

a  satellite  along  with  the  traraaitter  and  receiver 
positions  to  calculated  the  expected  Doppler  and 
chirp.  This  code  was  specifically  written  for  the 
bi-static  geometry  of  NAVSPASUR  at  216.98  MHz.  It 
can  be  easily  generalized  by  adding  FEMIT  to  the 
call. 

PROGRAM  ALGORITHM  (PSEUDOCODE):  The  actual  calculation  is  the  classic, 

non-relativistic  Doppler  shift.  The  calculation  is  repeated 
at  T=1  sec  to  obtain  the  chirp.  Acceleration  from  a  point 
source  Earth  and  corrections  for  rotating  coordinate  system 
are  included. 

INPUTS  EXPLICIT:  All  containing  in  call  list. 

OUTPUTS  EXPLICIT:  All  contained  in  call  list. 

MAJOR  VARIABLES: 

X,Y,Z,VX,VT,VZ  are  the  position  and  velocity  of  the  satellite  in  an 
Earth  centered  coordinate  system;  uiits  ■  and  n/sec. 

DOPLR, CHIRP  are  the  calculated  Doppler  and  chirp 

XR,YR,ZR,XT,YT,ZT  are  the  positions  of  the  reciver(R)  and  transmitter(T) 
in  Earth  centered,  rotating  coordinate  system. 

FEMIT  is  the  frequency  maittsd  by  the  transmitter.  For  NAVSPASUR, 

FEMIT  «  216.98  MHz. 

GM  is  gravitational  constant  for  the  Eerth. 

ONEGA  is  angular  velocity  of  the  Earth. 


MODIFIED: 


IMPLICIT  REAL'S  (A-Z) 

Enter  values  for  constants. 

FEMIT  >  216980000. 

VLCHT  ■  299792458. 

GM  ■  3.986124EU 

OMEGA  -  6.28318S3/(  23.*3600.  '  56.*60.  ) 

GM  is  the  gravitational  constant  and  OMEGA  is  the  angular  velocity  of 
the  eerth. 

Calculate  the  acceleration  of  the  satellite  based  on  position  and 
point  source  gravity  with  centrifigel  and  coriolia  force  added. 

Gravity  terms. 
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c 

GCNST  *  -I.O'GM/C  (  X*X  ♦  r*T  ♦  2*2  )  **1.5) 

AGX  X  GCNST  *  X 

AGY  *  GCNST  *  Y 

AG2  X  GCNST  *  2 

Coriolis  terms. 

ACORX  =  2*0HEGA*VY 
ACORY  X  -2*0MEGA*VX 

Centrifigal  terms. 

ACENX  X  OMEGA*QMEGA*X 
ACENY  X  QHEGA*ONEGA*Y 

Total  acceleration. 

AX  X  aGX  «  ACORX  +  ACENX 
AY  X  aGY  ♦  ACORY  ♦  ACENY 
A2  X  aG2 

Velocity  at  time  T=1  sec. 

VX1  X  vx  ♦  AX 

VY1  X  VY  ♦  AY 

V21  X  V2  ♦  A2 

Position  at  Txisec. 

XI  X  X  ♦  VX  ♦  .5*AX 

Y1  X  Y  ♦  VY  ♦  .5*AY 

21  X  2  ♦  V2  ♦  .5*A2 

Calculate  the  position  vector  and  distance  from  the  transmitter  to  the 

satellite  at  T*0  and  T«1sec. 

XSATx  X  •  XT 
YSATx  Y  -  YT 
2SAT»  2  •  2T 

XSAT1  X  XI  -  XT 
YSATI  X  Y1  .  YT 
2SAT1  X  Z1  -  2T 

D1  X  SORT<  XSAT*XSAT  ♦  YSAT*YSAT  ♦  2SAT*2SAT  ) 

05  X  S0RT<  XSAT1*XSAT1  ♦  YSAT1*TSAT1  ♦  2SAT1*2SAT1  ) 

Take  the  dot  product  of  the  velocity  and  the  unit  radius  vector  to 
get  velocity  along  line  of  sight. 

OOT1  -  (  XSAT*VX  ♦  ySAT*Vy  ♦  28AT*V2  )/01 
DOTS  X  <XSAT1*VX1  ♦  YSAT1*VY1  ♦  2SAT1*V21  )/05 

Calculate  the  Doppler  shifted  frequency. 

FREfll  X  FEhIT*<1.0  -  OOTlmCKT) 

FREQS  X  FENIT*(1.0  -  DOTS/VL6HT) 

Repeat  the  above  calculation  for  the  receiver. 

XSTA  X  X  •  XR 

YSTA  X  y  .  YR 

2STA  X  2  -  ZR 

XSTA1  X  XI  •  XR 

YSTA1  X  Y1  -  YR 

2$TA1  X  zl  -  2R 

C 

02  X  SQRTl  XSTA**2  ♦  YSTA**2  ♦  2STA**2  ) 
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D6  •  SQRT(  XSTA1*XSTA1  ♦  YSTA1*YSTA1  ♦  ZSTA1*ZSTA1) 

OOT2  »  <  XSTA*VX  ♦  YSTA*VY  ♦  ZSTA*VZ  )/D2 

D0T6  =  (XSTA1*VX1  ♦  YSTA1»VY1  ♦  2STA1*V21  )/D6 

FRE02  *  FREQ1*(100  -  00T2/VLCHT) 

FRE06  3  FREQS* (IDO  -  DOT6/VLGHT) 

Doppler  is  final  frequency  minus  transmit  frequency. 

**NOTE**  This  agrees  with  NAVSPASUR  sign  convention. 

Chirp  is  rate  of  change  of  Doppler. 

D(X>LR  s  FREQ2  •  FEMIT 
CHIRP  3  FRE06  -  FREa2 

RETURN 

END 


ooooooooonn  noono  ooriooooooooo 


SUBROUTINE  GAUSSJ(A.  N,  NP,  B,  N.  MP) 


This  routine  has  been  lifted  intact  froai  Nuaerical  Recipes. 

Consult  the  book  for  information  on  this  algorithm. 

LINEAR  EQUATION  SOLUTION  BY  GAUSS-JORDAN  ELIMINATION.  A  IS  AN 
INPUT  MATRIX  OF  N  BY  N  ELEMENTS,  STORED  IN  AN  ARRAY  OF  PHYSICAL 
DIMENSIONS  NP  BY  NP.  B  IS  AN  INPUT  MATRIX  Of  N  BY  M  CONTAINING 

THE  M  RIGHT-HAND  SIDE  VECTORS.  STORED  IN  AN  ARRAY  OF  PHYSICAL 
DIMENSIONS  NP  BY  NP.  ON  OUTPUT,  A  IS  REPLACED  BY  ITS  MATRIX 
INVERSE,  AND  B  IS  REPLACED  BY  THE  CORRESPONDING  SET  OF  SOLUTION 
VECTORS. 

IMPLICIT  REAL*8  (A-H.O-Z) 

PARAMETER(NMAX<50) 

DIMENSION  A(NP,NP).B(NP,MP),IPIV(NMAX),INDXR(NNAX),INOXC(NMAX) 

THE  INTEGER  ARRAYS  IPIV,  INOXR,  AND  INDXC  ARE  USED  FOR  BOOKKEEPING 
ON  THE  PIVOTING.  NMAX  SHOULD  BE  AS  LARGE  AS  THE  URGEST  ANTICI¬ 
PATED  VALUE  OF  N. 

DO  11  J>1.  N 
IPIV<J)*0 

11  CONTINUE 

DO  22  I>1.N 
BIG  >  0. 

DO  13  J>1.N 

IF(IPIV(J).NE.1)  THEN 
DO  12  K>1,N 

IF(IPIV(K).EQ.O)  THEN 

1F(ABS<A(J,K)).G£.81G)  THEN 
BIGaBS(A(J,K)) 

IROU-J 

ICOL>K 

END  IF 
ELSE 

IF(IPIV(K).GT.1}  STOP  'SINGULAR  MATRIX.' 

ENOIF 

12  CONTINUE 
ENDIF 

13  CONTINUE 

IPIV(ICOl)>IPIV(ICOL}  *■  1 

UE  NOU  HAVE  THE  PIVOT  ELEMENT,  SO  UE  INTERCHANGE  ROWS,  IF  NEEDED, 
TO  PUT  THE  PIVOT  ELEMENT  ON  THE  DIAGONAL.  THE  COLUMNS  ARE  NOT 
PHYSICALLY  INTERCHANGED,  ONLY  RELABELED:  INOX(I),  THE  COLUMN  OF 
THE  ITH  PIVOT  ELEMENT,  IS  THE  ITH  ELEMENT  THAT  IS  REDUCED,  WHILE 
INDXR(I)  IS  THE  ROW  IN  WHICH  THAT  PIVOT  ELEMENT  WAS  ORIGINALLY 
LOCATED.  IF  INOXR(I)  .HE.  INDXC(I)  THERE  IS  AN  IMPLIED  COLUMN 
INTERCHANGE.  WITH  THIS  FORM  OF  BOOKKEEPING,  THE  SOLUTION  B'S 
HILL  END  UP  IN  THE  CORRECT  ORDER,  AND  THE  INVERSE  MATRIX  HILL  BE 
SCRAMBLED  BY  COLUMNS. 

IF  (IROH.NE.ICOL)  THEN 
DO  U  L-1,H 

OUM«A(IROH,L) 

A(IROH,L)«A(ICOL.L> 

A(ICOL,L}>OUM 

U  CONTINUE 

DO  IS  L>1,M 

DUM>l<IROH,L) 

B(IROH,L)-B(ICOL,L) 

B(ICOL,L)«DUM 

15  CONTINUE 

ENOIF 

INOXR(I}*IROH 

INOXC(I)*ICOL 

IF(A(ICOL,ICOL).EQ.O.)  STOP  'SINGULAR  MATRIX' 

PIVINV  «  1./A(ICOL,ICOL) 

A(ICOL,ICOL)-1. 


SB 


16 


17 


18 


19 

21 

22 


23 

24 


DO  16  L>1.N 

A(ICOL.L)>A(ICOL.L>*PIVINV 
CONTINUE 
00  17  L«1,H 

B(ICOL,L)sB(ICOL.L)*PlVINV 

CONTINUE 

DO  21  LL>1,N 

IF(LL.NE.ICOL)  THEN 

OUN=A(LL.ICOL) 

A(LL,ICOL)*0. 

DO  18  L>1,H 

A(LL,L)«A<LL.L)-A(ICOL,L)*DUM 
CONTINUE 
DO  19  L>1.N 

B(LL.L>>B(LL.L)-B(IC0L,L)*DUN 

CONTINUE 

END  IF 
CONTINUE 

CONTINUE 

DO  24  L=N,1,-1 

IF(INOXR(L).NE.INDXC(L))  THEN 
DO  23  K>1,N 

OUM^AdC.INOXRCL)) 
A(K,INDXR(L))M(K.INOXC(L)) 
A(K.INOXC(t.))=OUH 
CONTINUE 
END  IF 
CONTINUE 

RETURN 

END 
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